library(lme4)
library(knitr)
library(survival)
library(coxme)
library(zoo)
library(ggthemes)
library(grid)
library(gridExtra)
library(ggbeeswarm)
library(stringr)
library(tidyverse)
library(parallel)
#Convenience functions
numextract <- function(string){
as.numeric(str_extract(string, "\\-*\\d+\\.*\\d*"))
}
comma <- function(x) format(x, digits = 3, big.mark = ",")
comma_1 <- function(x) format(x, digits = 1, big.mark = ",")
comma_2 <- function(x) format(x, digits = 2, big.mark = ",")
comma_p <- function(x){
if (x < 0.001){
return("< 0.001")
}
if (x<0.01 & x>=0.001){
paste("=", format(x, digits = 3, big.mark = ","))
}
else{
paste("=", format(x, digits = 2, big.mark = ","))
}
}
to_model <- read.csv("clean_time_series_revise2006_2015.csv")
#for table S4
tx_hr <- haven::read_sas("SAF 2018 Q3/tx_hr.sas7bdat", NULL) %>%
haven::zap_formats() %>% haven::zap_labels()
to_model
#three-status, donor quality, and transplant year
three_status <- "Surv(t_1, t_2, dead) ~ tx*(stat_2 + stat_1b)
+ tx_risky_don + era_tx"
#six-status, donor quality and transplant year
six_status <- "Surv(t_1, t_2, dead) ~ tx*(status_2 + status_3 + status_4 + status_6) +
tx_risky_don + era_tx"
formulas <- c(three_status,six_status)
fe_models <- vector(mode = "list", length(formulas))
hz_list <- vector(mode = "list", length(formulas))
for (i in seq_along(formulas)) {
coxph_fe <- coxph(as.formula(formulas[[i]]), to_model)
print(formulas[[i]])
print(coxph_fe)
fe_models[[i]] <- coxph_fe
#estimate baseline hazard function when all covariates equal to zero
# Status 1A, after "Status 1A", low risk donor
basehz <- basehaz(coxph_fe, centered=FALSE) %>%
mutate(hz_t = if_else(time ==1, hazard, (hazard - lag(hazard))/(time -lag(time))),
hz_t_time = if_else(time ==1, hz_t*time, hz_t*(time -lag(time))),
cum_hz = cumsum(hz_t_time)) %>% select(time, hz_t)
hz_list[[i]] <- basehz
}
me_models <- vector(mode = "list", length(formulas))
for (i in seq_along(formulas)) {
f_coxme <- as.formula(paste(formulas[[i]], " + (1+ tx|center)"))
start_time <- Sys.time()
three_status <- coxme(f_coxme, data = to_model)
me_models[[i]] <- three_status
end_time <- Sys.time()
diff_time <- end_time - start_time
print(f_coxme)
print(diff_time)
}
#save.image("model_fit.RData")
Here we write a function to calculate the between center combination in the survival benefit of transplant on the log hazard scale
\(survivalBenefit_i = \nu_{1i} - \nu_{0i}\)
\(Var(\nu_{i1} - \nu_{i0}) = Var(v_{1i}) + Var(\nu_{0i}) - 2*Cov(v_{1i}, v_{0i})\)
and finally since coxme returns the correlation between the random effects
\(Cov(v_{1i}, v_{0i}) = Corr(v_{1i}, v_{0i})*\sigma_{v_{0i}}\sigma_{v_{1i}}\)
variance_survival_benefit <- function(me_cox){
variance_est <- me_cox$vcoef[[1]]
var_wait <- variance_est[[1,1]]
var_tx <- variance_est[[2,2]]
cov_tx_wait <- variance_est[[1,2]]*sqrt(variance_est[[1,1]])*sqrt(variance_est[[2,2]])
return(var_wait + var_tx - 2*cov_tx_wait)
}
coxme resultsthree_status <- me_models[[1]]
basehz <- hz_list[[1]]
three_status_var <- variance_survival_benefit(three_status)
three_status
Cox mixed-effects model fit by maximum likelihood
Data: to_model
events, n = 11058, 163240
Iterations= 30 184
NULL Integrated Fitted
Log-likelihood -107600 -105890.7 -105702
Chisq df p AIC BIC
Integrated loglik 3418.65 10.00 0 3398.65 3325.54
Penalized loglik 3795.97 122.76 0 3550.45 2652.96
Model: f_coxme
Fixed coefficients
coef exp(coef) se(coef) z p
tx -1.95316608 0.1418243 0.04469987 -43.70 0.000
stat_2 -1.35517243 0.2579028 0.03648152 -37.15 0.000
stat_1b -0.94242273 0.3896826 0.03167335 -29.75 0.000
tx_risky_don 0.26049398 1.2975709 0.02801652 9.30 0.000
era_tx -0.05394094 0.9474881 0.03047029 -1.77 0.077
tx:stat_2 1.32075844 3.7462616 0.06369779 20.73 0.000
tx:stat_1b 0.90417441 2.4698920 0.04370071 20.69 0.000
Random effects
Group Variable Std Dev Variance Corr
center Intercept 0.22364632 0.05001768 -0.50908712
tx 0.23822093 0.05674921
The calculated variance in the survival benefit of transplant \(B = Tx - Wait\) is 0.161 on log hazard ratio scale for the 3-status model
coxme resultssix_status <- me_models[[2]]
basehz_six_stat <- hz_list[[2]]
six_status_var <- variance_survival_benefit(six_status)
six_status
Cox mixed-effects model fit by maximum likelihood
Data: to_model
events, n = 11058, 163240
Iterations= 23 211
NULL Integrated Fitted
Log-likelihood -107600 -105838.3 -105664.7
Chisq df p AIC BIC
Integrated loglik 3523.39 14.00 0 3495.39 3393.04
Penalized loglik 3870.57 124.95 0 3620.67 2707.17
Model: f_coxme
Fixed coefficients
coef exp(coef) se(coef) z p
tx -2.92406175 0.05371507 0.13932787 -20.99 0.0e+00
status_2 -0.86879490 0.41945673 0.08926304 -9.73 0.0e+00
status_3 -1.52574756 0.21745843 0.08341592 -18.29 0.0e+00
status_4 -2.16907702 0.11428305 0.07852884 -27.62 0.0e+00
status_6 -2.61411255 0.07323275 0.08224010 -31.79 0.0e+00
tx_risky_don 0.26222497 1.29981894 0.02800619 9.36 0.0e+00
era_tx -0.05269387 0.94867038 0.03042904 -1.73 8.3e-02
tx:status_2 0.73718240 2.09003832 0.15245228 4.84 1.3e-06
tx:status_3 1.24972334 3.48937744 0.14226535 8.78 0.0e+00
tx:status_4 1.85129304 6.36804832 0.13835133 13.38 0.0e+00
tx:status_6 2.28656875 9.84111238 0.14941169 15.30 0.0e+00
Random effects
Group Variable Std Dev Variance Corr
center Intercept 0.19542629 0.03819143 -0.38557673
tx 0.21220860 0.04503249
The calculated variance in the survival benefit of transplant \(B = Tx - Wait\) in the six-status model decreased by a factor of 0.28 from 0.161 to 0.115 on log hazard ratio scale compared to the three-status system.
hr_tx_ij <- function(model = three_status, df= to_model, i, j){
#this function returns the log hazard ratio of transplant
#for patient j transplanted at center i for a given coxme model
fixed_effects <- fixef(model)
x_vars <- names(fixed_effects)
ref <- data.frame(ranef(model)[[1]])
ref$center <- row.names(ref)
tx_ref <- filter(ref, center == i)$tx
is_int <- function(x) grepl("tx:", x)
ints <- x_vars[sapply(x_vars, is_int)]
non_ints <- setdiff(x_vars, ints)
d_ij <- df %>% filter(PX_ID == j) %>% group_by(tx) %>%
filter(row_number()==1 & tx ==1) %>% select(one_of(non_ints))
h_ij <- tx_ref + unname(fixed_effects['tx'])
for (x in ints) {
fe_tx <- unname(fixed_effects[x])
c_list <- strsplit(x, ":")
int_var <- c_list[[1]][2]
h_ij <- h_ij + d_ij[[int_var]]*fe_tx
}
return(h_ij)
}
Warnings:
b_ij <- function(df, base_hz, model, i, j){
fixed_effects <- fixef(model)
x_vars <- names(fixed_effects)
ref <- data.frame(ranef(model)[[1]])
ref$center <- row.names(ref)
d_ij <- df %>% filter(PX_ID == j) %>%
select(PX_ID, t_1, t_2, dead, one_of(x_vars)) %>%
mutate(time = t_1,
time = ifelse(time==0, 1, time),
fe_lp = 0, lp_no_tx = 0)
for (x in x_vars) {
fe_x <- unname(fixed_effects[x])
if(grepl("tx:", x) == FALSE){
d_ij$c_x <- d_ij[[x]]
}
##add code to deal with interaction terms here
if(grepl("tx:", x) == TRUE){
c_list <- strsplit(x, ":")
int_var <- c_list[[1]][2]
d_ij$c_x <- d_ij[[int_var]]*d_ij[["tx"]]
}
d_ij <- d_ij %>% mutate(fe_lp = fe_lp + c_x*fe_x) %>% select(-c_x)
}
###Hypothetical when patient has always been transplanted
non_tx_terms <- lapply(x_vars, function(x) if(grepl("tx", x) == FALSE) x)
non_tx_terms <- non_tx_terms[-which(sapply(non_tx_terms, is.null))]
for (x in non_tx_terms) {
fe_x <- unname(fixed_effects[x])
d_ij$c_x <- d_ij[[x]]
d_ij <- d_ij %>% mutate(lp_no_tx = lp_no_tx + c_x*fe_x) %>% select(-c_x)
}
d_ij <- left_join(base_hz, d_ij, by = c("time")) %>%
select(-t_1, -t_2)
d_ij <- d_ij %>%
mutate_all(funs(na.locf(., na.rm = FALSE)))
ctr_i_b0 <- ref[ref$center == i,]$Intercept
ctr_i_b1 <- ref[ref$center == i,]$tx
d_ij <- d_ij %>% mutate(
fe_wait = lp_no_tx,
fe_tx = fe_lp,
ctr_hz_wait = hz_t*exp(fe_wait + ctr_i_b0),
ctr_hz_tx = hz_t*exp(fe_tx + ctr_i_b0 + ctr_i_b1)
) %>%
mutate(wait_surv = cumprod(1-ctr_hz_wait)) %>%
group_by(tx) %>%
mutate(
surv_post_tx = if_else( tx ==1, cumprod(1-ctr_hz_tx), NULL),
surv_post_tx = max(wait_surv)*surv_post_tx) %>% ungroup()
tx_point <- d_ij %>%
group_by(tx) %>%
filter(row_number()==1 & tx==1) %>%
ungroup()
if (nrow(tx_point)==0) warning(paste0("Candidate ", j, " not transplanted"))
if (nrow(tx_point)>0) {
plot_i_j <- ggplot() +
geom_step(data = d_ij,
aes(y =wait_surv,
x = time,
color = "Waitlist Survival"),
size = 1.5) +
geom_step(data = d_ij,
aes(y =surv_post_tx,
x = time,
color = "Post-transplant survival"),
size = 1.5) +
geom_ribbon(data = d_ij,
aes(x= time,
ymin = wait_surv,
ymax = surv_post_tx,
fill = "Survival Benefit"),
alpha = 0.5) +
geom_point(data = tx_point,
aes(x = time,
y = surv_post_tx,
shape= "Transplant"),
size = 5, colour = "darkgreen") +
labs(x = "Time(days)", y = "Survival", shape = "", fill = "") +
scale_y_continuous(limits = c(0,1), breaks = seq(0,1, 0.25)) +
scale_x_continuous(breaks = seq(0, 3650, 365)) +
scale_color_manual(name = NULL, values=c("blue", "red")) +
scale_fill_manual(name = "", values = c("lightblue"))
return(plot_i_j)
}
}
this function calculates the survival benefit for a given set of times (default is five-years). Calculates the survival function conditional upon reaching transplant (“resets” the survival curve to 100% at transplant)
Warnings:
ben_wait_tx <- function(df = to_model,
base_hz = basehz,
model = three_status, i, j, end_times = 1825){
max_end_time <- max(end_times)
fixed_effects <- fixef(model)
x_vars <- names(fixed_effects)
ref <- data.frame(ranef(model)[[1]])
ref$center <- row.names(ref)
d_ij <- df %>%
filter(PX_ID == j) %>%
select(PX_ID, t_1, t_2, dead, one_of(x_vars)) %>%
mutate(time = t_1,
time = ifelse(time==0, 1, time),
fe_lp = 0,
lp_no_tx = 0)
for (x in x_vars) {
fe_x <- unname(fixed_effects[x])
if(grepl("tx:", x) == FALSE){
d_ij$c_x <- d_ij[[x]]
}
##add code to deal with interaction terms here
if(grepl("tx:", x) == TRUE){
c_list <- strsplit(x, ":")
int_var <- c_list[[1]][2]
d_ij$c_x <- d_ij[[int_var]]*d_ij[["tx"]]
}
d_ij <- d_ij %>% mutate(fe_lp = fe_lp + c_x*fe_x) %>% select(-c_x)
}
###Hypothetical when patient has always been transplanted
non_tx_terms <- lapply(x_vars, function(x) if(grepl("tx", x) == FALSE) x)
non_tx_terms <- non_tx_terms[-which(sapply(non_tx_terms, is.null))]
for (x in non_tx_terms) {
fe_x <- unname(fixed_effects[x])
d_ij$c_x <- d_ij[[x]]
d_ij <- d_ij %>% mutate(lp_no_tx = lp_no_tx + c_x*fe_x) %>% select(-c_x)
}
d_ij <- left_join(base_hz, d_ij, by = c("time")) %>% select(-t_1, -t_2)
d_ij <- d_ij %>%
mutate_all(funs(na.locf(., na.rm = FALSE))) %>%
filter(tx ==1) %>%
mutate(time = time - min(time)+1)
ctr_i_b0 <- ref[ref$center == i,]$Intercept
ctr_i_b1 <- ref[ref$center == i,]$tx
d_ij <- d_ij %>% mutate(
fe_wait = lp_no_tx,
fe_tx = fe_lp,
ctr_hz_wait = hz_t*exp(fe_wait + ctr_i_b0),
ctr_hz_tx = hz_t*exp(fe_tx + ctr_i_b0 + ctr_i_b1)
) %>%
mutate(wait_surv = cumprod(1-ctr_hz_wait)) %>%
group_by(tx) %>%
mutate(
surv_post_tx = if_else( tx ==1, cumprod(1-ctr_hz_tx), NULL),
surv_post_tx = max(wait_surv)*surv_post_tx) %>% ungroup()
d_ij <- d_ij %>%
filter(time <= max_end_time) %>%
mutate(surv_benefit = surv_post_tx - wait_surv)
if (nrow(d_ij)==0) warning(paste0("Candidate ", j, " not transplanted"))
s_benefit <- vector(mode = "numeric", length = 2)
wait_5 <- filter(d_ij, time == end_times)$wait_surv
post_5 <- filter(d_ij, time == end_times)$surv_post_tx
if (length(wait_5) == 0){
wait_5 <- NA
}
if (length(post_5)==0){
post_5 <- NA
}
s_benefit[[1]] <- wait_5
s_benefit[[2]] <- post_5
return(s_benefit)
}
Warnings:
detectCores()
no_cores <- detectCores() - 2
# Initiate cluster
cl <- makeCluster(no_cores, type="FORK")
center_loop <- function(centers, recips, cl){
#output tibbles
wait_surv <- tibble(recips = recips)
post_tx_surv <- tibble(recips = recips)
for (c in centers) {
cur_center <- as.numeric(c)
print(paste0("currently calculating for center ", c))
B_i <- parSapply(cl, r_list,
function(x) ben_wait_tx(j = x,
df = to_model,
base_hz = basehz,
model = three_status,
i = cur_center))
B_trans <- matrix(B_i, nrow = 2, ncol = length(recips)) %>% t()
B_tib <- as.tibble(B_trans)
wait <- B_tib %>% select(V1) %>% cbind(recips)
colnames(wait) <- c(c, "recips")
post_tx <- B_tib %>% select(V2) %>% cbind(recips)
colnames(post_tx) <- c(c, "recips")
wait_surv <- wait_surv %>% left_join(wait, by = "recips")
post_tx_surv <- post_tx_surv %>% left_join(post_tx, by = "recips")
}
benefit <- list(wait_surv, post_tx_surv)
return(benefit)
}
r_list <- unique(filter(to_model, tx ==1)$PX_ID)
c_list <- as.numeric(row.names(ranef(three_status)[[1]]))
by_time_tibbles <- center_loop(c_list, r_list, cl)
stopCluster(cl)
five_yr_wait <- by_time_tibbles[[1]]
five_yr_post <- by_time_tibbles[[2]]
#recipients
recips <- to_model %>% group_by(PX_ID) %>% filter(tx == 1) %>% ungroup()
#marginal distribution of benefit of transplant
marginal_hazards <- mapply(hr_tx_ij, i = recips$center, j = recips$PX_ID)
marginal_waits_posts <- mapply(ben_wait_tx, i = recips$center, j = recips$PX_ID)
marginal_waits <- marginal_waits_posts[1,]
marginal_posts <- marginal_waits_posts[2,]
pop_tx_hz <- recips %>% cbind(marginal_hazards, marginal_waits, marginal_posts) %>%
mutate(log_hazard = marginal_hazards,
hr = exp(log_hazard),
benefit= marginal_posts- marginal_waits)
r_list <- unique(filter(to_model, tx ==1)$PX_ID)
c_list <- as.numeric(row.names(ranef(three_status)[[1]]))
center_OPO <- recips %>% select(center, OPO) %>%
group_by(center, OPO) %>%
summarise(num_tx = n()) %>% ungroup()
re_effects <- data.frame(ranef(three_status)[[1]])
re_effects$center <- row.names(re_effects)
colnames(re_effects)[2] <- "tx_re"
#make a center level dataset of random effects
center_re <- center_OPO %>%
left_join(re_effects %>% mutate(center = as.numeric(center))) %>%
mutate(tx_exp = exp(tx_re),
benefit = tx_re-Intercept,
intercept_exp = exp(Intercept)) %>% ungroup()
center_re <- center_re %>% left_join(five_year_means)
#calculate the minimum and maximum center
min_tx_exp <- min(center_re$tx_exp)
max_tx_exp <- max(center_re$tx_exp)
min_int_exp <- min(center_re$intercept_exp)
max_int_exp <- max(center_re$intercept_exp)
#make benefit dataframe- hypothetical estimates as if all recipients were transplanted at the center
colnames(five_yr_wait) <- c("recips", c_list)
colnames(five_yr_post) <- c("recips", c_list)
M <- merge(five_yr_post,five_yr_wait,by="recips")
ben_five_year_only <- M[,grepl("*\\.x$",names(M))] - M[,grepl("*\\.y$",names(M))]
cbind(M[,1,drop=FALSE],ben_five_year_only)
remove(M)
five_year_waits <- five_yr_wait %>% select(-recips) %>%
summarise_all(mean, na.rm =TRUE) %>%
gather(value = mean_five_year_wait, key = center) %>%
select(center, mean_five_year_wait) %>% mutate(center = as.numeric(center))
five_year_posts <- five_yr_post %>% select(-recips) %>%
summarise_all(mean, na.rm =TRUE) %>%
gather(value = mean_five_year_post, key = center) %>%
select(center, mean_five_year_post) %>% mutate(center = as.numeric(center))
#center level dataset of summary statistics for case-mix adjusted
#five-year mean waitlist, post-transplant, and survival benefits
five_year_means <- ben_five_year_only %>%
summarise_all(funs(mean,
median,
sd,
IQR,
lowQR = quantile(., probs=0.25),
upQR = quantile(., probs=0.75)), na.rm = TRUE) %>%
gather(key = "key", value = "value") %>%
separate(key, c("center", "stat"), sep = "_") %>%
spread(stat, value) %>%
mutate(center = numextract(center), mean_five_year = mean) %>% select(-mean) %>%
left_join(five_year_waits, by = "center") %>%
left_join(five_year_posts, by = "center") %>%
left_join(center_OPO, by = "center")
five_year_means <- five_year_means %>%
mutate(se_mean = sd/sqrt(num_tx),
upper_ci = mean_five_year + 1.96*se_mean,
lower_ci = mean_five_year - 1.96*se_mean,
sig_better = if_else(lower_ci >mean(five_year_means$mean_five_year), 1, 0),
sig_worse = if_else(upper_ci <mean(five_year_means$mean_five_year), 1, 0),
center_performance = case_when(
lower_ci >mean(five_year_means$mean_five_year) ~ "High",
upper_ci <mean(five_year_means$mean_five_year) ~ "Low",
TRUE ~ "Average"
), center_performance = factor(center_performance,
levels = c("Low", "Average", "High")))
#marginal distribution of benefit of transplant under six status model
marginal_hazards_6_stat <- mapply(hr_tx_ij,
i = recips$center,
j = recips$PX_ID,
MoreArgs = list(model = six_status))
marginal_waits_posts_6_stat <- mapply(ben_wait_tx,
i = recips$center,
j = recips$PX_ID,
MoreArgs = list(base_hz = basehz_six_stat,
model = six_status))
marginal_waits_6_stat <- marginal_waits_posts_6_stat[1,]
marginal_posts_6_stat <- marginal_waits_posts_6_stat[2,]
pop_tx_hz_6_stat <- recips %>% cbind(marginal_waits_6_stat, marginal_posts_6_stat) %>%
mutate(benefit_6 = marginal_posts_6_stat - marginal_waits_6_stat)
#not run- save workspace image
#save.image("model_fit_06_15_w_benefits.RData")
to_plot <- five_year_means %>% arrange(mean_five_year) %>% mutate(index = row_number())
g_legend<-function(a.gplot){
tmp <- ggplot_gtable(ggplot_build(a.gplot))
leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box")
legend <- tmp$grobs[[leg]]
return(legend)}
catepillar <- ggplot(to_plot, aes(x = index,
y = mean_five_year*100,
ymin = lower_ci*100,
ymax = upper_ci*100,
color = center_performance)) +
geom_point() +
geom_errorbar() +
scale_color_manual(name = "Center benefit:",
values= c("red3","black","navyblue"),
guide = guide_legend()) +
scale_x_continuous(breaks = c(1, 113)) +
scale_y_continuous(breaks = seq(25, 60, 5), limits = c(25, 60)) +
labs(y = "Absolute five-year survival benefit \nof heart transplant (%)",
x = "Center", color = "") +
theme_few() +
theme(axis.title = element_text(),
axis.title.y = element_text(),
legend.title=element_text(size=12),
legend.text=element_text(size=12),
panel.grid.major.y = element_line(color = "gray", size = 0.25),
legend.position="bottom")
box <- ggplot(to_plot %>% mutate(all = " "), aes(x = all, y = 100*mean_five_year)) +
geom_boxplot() +
stat_boxplot(geom = "errorbar", width = 0.5) +
geom_point(y = mean(to_plot$mean_five_year)*100,
aes(x =all, shape = "mean"),
size=4, fill = "gray", alpha = 0.5 )+
theme_few() +
labs(x = " ", shape = " ") +
scale_shape_manual(values = 23) +
theme(
axis.title.y = element_blank(),
axis.ticks.y = element_blank(),
axis.text.y = element_blank(),
legend.title=element_text(size=12),
legend.text=element_text(size=12),
panel.grid.major.y = element_line(color = "gray", size = 0.25),
legend.position= "bottom",
legend.justification='left') +
scale_y_continuous(limits = c(25, 60), breaks = seq(25, 60, 5))
fig_2 <- grid.arrange(catepillar, box,
nrow=1, widths = c(1, 0.25))
num_high <- table(five_year_means$center_performance)["High"][[1]]
num_low <- table(five_year_means$center_performance)["Low"][[1]]
max_benefit <- comma_2(100*max(to_plot$mean_five_year))
min_benefit <- comma_2(100*min(to_plot$mean_five_year))
high_qrt_5 <- comma_2(100*quantile(to_plot$mean_five_year, probs = 0.75))
low_qrt_5 <- comma_2(100*quantile(to_plot$mean_five_year, probs = 0.25))
median_5 <- comma_2(100*median(to_plot$mean_five_year))
ggsave( "fig_2.pdf", plot = fig_2, width = 6, height = 4)
Between-center variation in the absolute improvement in mean five-year survival associated with heart transplantation. Caterpillar plot (left) of center estimates case-mix adjusted for time from listing to transplant, waitlist Status, year of listing, donor quality, and cold ischemic time. Ninety-five percent confidence intervals are constructed using the variance of the distribution of five-year survival estimates for the entire study population and the number of transplants performed at each center during the study period. Compared to the mean center, 31 centers (27%) had significantly higher survival benefit and 30 centers (27%) had lower benefit.
The adjusted five-year survival benefit varied substantially by center, ranging from 30% to 55%, median 43 and IQR [40% - 47%] (Box plot, right). Whiskers represent smallest and largest observations within 1.5 times the IQR of the hinges, points represent outliers beyond this range. The mean of center five-year survival benefits is represented by a gray diamond.
panel_A <- ggplot(data = center_re,
aes(x = 100*mean_five_year_wait,
y = 100*mean_five_year)) +
geom_point(aes(color = center_performance), size = 2.5, alpha = 0.9) +
#annotate("label", x = 25, y = 30, label = paste0("r = ", cor_model), size = 5) +
scale_x_continuous(name = "Five-year survival without transplant (%)",
limits = c(18, 44),
breaks = c(20, 25, 30, 35, 40)) +
scale_y_continuous(name = "Center five-year surival benefit (%)",
limits = c(27, 60), breaks = seq(30, 60, 5)) +
scale_color_manual(name = "Center benefit:",
values = c("red3","black","navyblue")) +
theme_few() +
ggtitle("Panel A") +
theme(legend.key = element_rect(fill = "white"),
legend.background = element_rect(fill = "white"),
legend.direction = "horizontal",
panel.border = element_blank(),
axis.line = element_line(color = "black", size = 0.25),
plot.title = element_text(hjust = 0.5),
legend.title=element_text(size=12),
legend.text=element_text(size=12),
legend.box.background = element_rect(colour = "black"))
mylegend <- g_legend(panel_A)
panel_B <- ggplot(data = center_re,
aes(x = mean_five_year_post*100, y = 100*mean_five_year)) +
geom_point(aes(color = center_performance), size = 2.5, alpha = 0.9) +
scale_x_continuous(name = "Five-year post-transplant survival (%)",
limits = c(63, 89),
breaks = c(65, 70, 75, 80, 85)) +
scale_y_continuous(limits = c(27, 60), breaks = seq(30, 60, 5)) +
scale_color_manual(values = c("red3","black","navyblue")) +
theme_few() +
ggtitle("Panel B")+
theme(axis.title.y = element_blank(),
panel.border = element_blank(),
axis.line = element_line(color = "black", size = 0.25),
legend.position="none",
panel.grid.major.y = element_line(color = "gray", size = 0.25),
plot.title = element_text(hjust = 0.5))
fig_3 <- grid.arrange(
arrangeGrob(panel_A + theme(legend.position="none",
panel.grid.major.y = element_line(color = "gray", size = 0.25)),
panel_B,
nrow=1, widths = c(1, 0.95)),
mylegend,
nrow=2,heights=c(10, 2.2))
slope_waitlist <- lm(mean_five_year ~ mean_five_year_wait,
data = center_re)
slope_post <- lm(mean_five_year ~ mean_five_year_post,
data = center_re)
pct_benefit_increase_per_neg_10 <- comma_2(
10*-slope_waitlist$coefficients[[2]])
ci_slope_wait_low <- comma_2(10*-confint(slope_waitlist, 'mean_five_year_wait', level=0.95)[1,1])
ci_slope_wait_up <- comma_2(10*-confint(slope_waitlist, 'mean_five_year_wait', level=0.95)[1,2])
pct_benefit_increase_post <- comma_2(10*slope_post$coefficients[[2]])
ci_slope_post_low <- comma_2(10*-confint(slope_post, 'mean_five_year_post', level=0.95)[1,1])
ci_slope_post_up <- comma_2(10*-confint(slope_post, 'mean_five_year_post', level=0.95)[1,2])
ggsave("fig_3.pdf", plot = fig_3)
Mean five-year survival without transplant, five-year post-transplant survival, and five-year survival benefit for each US heart transplant center. There was a significant relationship between the Status-adjusted medical urgency of candidates listed by each center (as measured by risk of death on the waitlist) and the benefit of transplantation measured by five-year survival benefit (left panel). For every 10% decrease in expected candidate waitlist surival, there was a 6.2% increase in estimated survival benefit associated with heart transplant (95% CI 5.2% - 7.3%). In contrast, there was no significant relationship between post-transplant survival and benefit (+1.5%, 95% CI -3.8%- 0.83%) (right panel).
five_yr_by_status <- pop_tx_hz %>%
group_by(three_status) %>%
summarise(
number = n(),
median = 100*median(benefit, na.rm = TRUE),
upper_IQR = 100*quantile(benefit, probs = 0.75, na.rm = TRUE),
lower_IQR = 100*quantile(benefit, probs = 0.25, na.rm = TRUE),
mean = 100*mean(benefit, na.rm = TRUE),
mean_waitlist = 100*mean(marginal_waits, na.rm = TRUE),
mean_post = 100*mean(marginal_posts, na.rm = TRUE)
)
five_yr_by_status
three_box <- ggplot(data = pop_tx_hz %>%
mutate(x_labels = case_when(
three_status == "Status 1A" ~ paste0("Status 1A\nN =", comma_2(filter(five_yr_by_status, three_status == "Status 1A")$number)),
three_status == "Status 1B" ~ paste0("Status 1B\nN =", comma_2(filter(five_yr_by_status, three_status == "Status 1B")$number)),
three_status == "Status 2" ~ paste0("Status 2\nN =", comma_2(filter(five_yr_by_status, three_status == "Status 2")$number))
)),
aes(y = 100*benefit, x=x_labels, color = three_status, fill = three_status)) +
geom_dotplot(binaxis = "y",
stackdir='center',
binwidth = 1, dotsize = 0.2, stackratio = .5) +
guides(color=FALSE, fill = FALSE) +
labs(y = "Absolute 5-year survival benefit (%)", fill = "") +
scale_y_continuous(limits = c(-5, 100)) +
theme_few() +
theme(axis.title.x = element_blank(),
panel.grid.major.y = element_line(color = "gray", size = 0.25)) +
scale_x_discrete(expand = expand_scale( mult = c(0.5, 0.2))) +
scale_color_brewer(palette = "Set1", direction = -1) +
scale_fill_brewer(palette = "Set1", direction = -1) +
stat_boxplot(aes(x = x_labels, y = 100*benefit),
geom = "errorbar",
color = "black",
width = 0.25, lwd=0.75) +
geom_boxplot(aes(x = x_labels, y = 100*benefit),
alpha = 0, color = "black",
outlier.shape = NA,
width = 0.5, lwd=0.75, fatten = 1.5)
three_box
ggsave("fig_4A.eps", plot = three_box, width = 8, height = 5)
Dot and box plots of estimated five-year survival benefit associated with heart transplantation for adult recipients by Status at transplantation between 2006-2015. There were 11,227 Status 1A recipients with median benefit from transplant of 59%, [IQR 54- 62%]. There were 7,250 Status 1B recipients with median benefit from transplant of 27%, [IQR 23- 31%]. In addition, there were only 1,338 Status 2 recipients with median benefit from transplant of 14%, [IQR 10- 17%]. Within priority Status groups, the standard deviation in the five-year survival benefit associated with transplant was 5.5% with the majority (76%) of the variation amongst recipients attributable to center-level effects and the remaining 24% attributable to variation in transplant year, waiting time, donor quality and ischemic time.
line_thickness <- 0.55
fatten_level <- 1.5
five_yr_by_six_status <- pop_tx_hz_6_stat %>%
group_by(six_status) %>%
summarise(
number = n(),
median = 100*median(benefit_6, na.rm = TRUE),
upper_IQR = 100*quantile(benefit_6, probs = 0.75, na.rm = TRUE),
lower_IQR = 100*quantile(benefit_6, probs = 0.25, na.rm = TRUE),
mean = 100*mean(benefit_6, na.rm = TRUE),
mean_waitlist = 100*mean(marginal_waits_6_stat, na.rm = TRUE),
mean_post = 100*mean(marginal_posts_6_stat, na.rm = TRUE)
)
g1 <- ggplot(pop_tx_hz_6_stat %>%
filter(six_status == 1) %>%
mutate(six_status = paste0("Status 1\nN =", comma_2(filter(five_yr_by_six_status, six_status == 1)$number)))) +
geom_dotplot(aes(x = six_status, y = 100*benefit_6),
color = RColorBrewer::brewer.pal(3, "Set1")[[3]],
fill = RColorBrewer::brewer.pal(3, "Set1")[[3]],
binaxis = "y", stackdir = "center",
method="dotdensity", stackgroups = T, binpositions="all",
binwidth = 1, dotsize = 0.2, stackratio = .5) +
theme_few() +theme(axis.title.x = element_blank(),
panel.border=element_blank(),
axis.line = element_line(color="black", size = 0.25),
legend.position = "none",
plot.margin = margin(0,-4,0,2),
panel.grid.major.y = element_line(color = "gray", size = 0.25)) +
scale_y_continuous(name = "Absolute 5-year survival benefit (%)", limits = c(-5, 100)) +
stat_boxplot(aes(x = six_status, y = 100*benefit_6),
geom = "errorbar", color = "black",
width = 0.25, lwd=line_thickness) +
geom_boxplot(aes(x = six_status, y = 100*benefit_6),
fill = "gray", alpha = 0, color = "black",
outlier.shape = NA, width = 0.5,
lwd=line_thickness, fatten = fatten_level)
for (i in c(2,3)){
assign(paste0("g", i), ggplot(pop_tx_hz_6_stat %>%
filter(six_status == i) %>%
mutate(six_status = case_when(
six_status == 1 ~paste0("Status 1\nN =", comma_2(filter(five_yr_by_six_status, six_status == 1)$number)),
six_status == 2 ~paste0("Status 2\nN =", comma_2(filter(five_yr_by_six_status, six_status == 2)$number)),
six_status == 3 ~paste0("Status 3\nN =", comma_2(filter(five_yr_by_six_status, six_status == 3)$number)),
six_status == 4 ~paste0("Status 4\nN =", comma_2(filter(five_yr_by_six_status, six_status == 4)$number)),
six_status == 6 ~paste0("Status 6\nN =", comma_2(filter(five_yr_by_six_status, six_status == 6)$number))
))) +
geom_dotplot(aes(x = six_status, y = 100*benefit_6),
color = RColorBrewer::brewer.pal(3, "Set1")[[3]], fill = RColorBrewer::brewer.pal(3, "Set1")[[3]],
binaxis = "y", stackdir = "center", method="dotdensity", stackgroups = T, binpositions="all",
binwidth = 1, dotsize = 0.2, stackratio = .5) +
labs(fill = "", color = "") +
theme_few() +theme(axis.title.x = element_blank(),
panel.border=element_blank(),
legend.position = "bottom",
axis.line = element_line(color="black", size = 0.25),
axis.line.y = element_blank(),
axis.ticks.y = element_blank(),
axis.title.y = element_blank(),
axis.text.y =element_blank(),
plot.margin = margin(0,-4,0,-4),
panel.grid.major.y = element_line(color = "gray", size = 0.25))+
scale_y_continuous(limits = c(-5, 100)) +
scale_color_brewer(palette = "Set1", direction = 1) +
scale_fill_brewer(palette = "Set1", direction = 1) +
stat_boxplot(aes(x = six_status, y = 100*benefit_6),
geom = "errorbar", color = "black",
width = 0.25, lwd=line_thickness) +
geom_boxplot(aes(x = six_status, y = 100*benefit_6),
fill = "gray", alpha = 0, color = "black",
lwd=line_thickness, fatten = fatten_level)
)
}
for (i in c(4, 6)){
assign(paste0("g", i), ggplot(pop_tx_hz_6_stat %>%
filter(six_status == i) %>%
mutate(six_status = case_when(
six_status == 1 ~paste0("Status 1\nN =", comma_2(filter(five_yr_by_six_status, six_status == 1)$number)),
six_status == 2 ~paste0("Status 2\nN =", comma_2(filter(five_yr_by_six_status, six_status == 2)$number)),
six_status == 3 ~paste0("Status 3\nN =", comma_2(filter(five_yr_by_six_status, six_status == 3)$number)),
six_status == 4 ~paste0("Status 4\nN =", comma_2(filter(five_yr_by_six_status, six_status == 4)$number)),
six_status == 6 ~paste0("Status 6\nN =", comma_2(filter(five_yr_by_six_status, six_status == 6)$number))
))) +
geom_dotplot(aes(x = six_status, y = 100*benefit_6, color = three_status, fill=three_status),
binaxis = "y", stackdir = "center",
method="dotdensity", stackgroups = T, binpositions="all",
binwidth = 1, dotsize = 0.2, stackratio = .5) +
labs(fill = "Status in 3 tier system:", color = "Status in 3 tier system:") +
theme_few() +theme(axis.title.x = element_blank(),
panel.border=element_blank(),
legend.position = "bottom",
axis.line = element_line(color="black", size = 0.25),
axis.line.y = element_blank(),
axis.ticks.y = element_blank(),
axis.title.y = element_blank(),
axis.text.y =element_blank(),
legend.title=element_text(size=11),
legend.text=element_text(size=11),
legend.background = element_blank(),
legend.box.background = element_rect(colour = "black"),
plot.margin = margin(0,-4,0,-4),
panel.grid.major.y = element_line(color = "gray", size = 0.25))+
scale_y_continuous(limits = c(-5, 100)) +
scale_color_brewer(palette = "Set1", direction = -1) +
scale_fill_brewer(palette = "Set1", direction = -1) +
stat_boxplot(aes(x = six_status, y = 100*benefit_6),
geom = "errorbar", color = "black",
width = 0.25, lwd=line_thickness) +
geom_boxplot(aes(x = six_status, y = 100*benefit_6),
fill = "gray", alpha = 0,
color = "black", outlier.shape = NA,
width = 0.5, lwd=line_thickness, fatten = fatten_level)
)
}
mylegend <- g_legend(g4)
Fig_4B <- grid.arrange(
arrangeGrob(g1 + theme(legend.position="none",
panel.grid.major.y = element_line(color = "gray", size = 0.25)),
g2 + theme(legend.position="none",
panel.grid.major.y = element_line(color = "gray", size = 0.25)),
g3 + theme(legend.position="none",
panel.grid.major.y = element_line(color = "gray", size = 0.25)),
g4 + theme(legend.position="none",
panel.grid.major.y = element_line(color = "gray", size = 0.25)),
g6 + theme(legend.position="none",
panel.grid.major.y = element_line(color = "gray", size = 0.25)),
nrow=1, widths = c(0.5, 0.3, 0.7, 1.2, 0.8)),
mylegend,
nrow=2,heights=c(10, 2.2))
ggsave("fig_4B.eps", plot = Fig_4B, width = 10, height = 5)
Dot and box plot of estimated five-year survival benefit associated with heart transplantation for adult recipients by projected six-status at transplantation between 2006-2015. Dot color corresponds to three-status at transplant, x-axis corresponds to re-assigned six-tier Status. There would have been 255 Status 1 recipients with median benefit from transplant of 69%, [IQR 65- 71%], 1,254 Status 2 recipients with median benefit from transplant of 65%, [IQR 61- 67%], 6,255 Status 3 recipients with median benefit from transplant of 49%, [IQR 46- 53%], 11,000 Status 4 recipients with median benefit from transplant of 28%, [IQR 24- 31%], and 1,051 Status 6 recipients with median benefit from transplant of 14%, [IQR 10- 17%]. Whiskers represent smallest and largest observations within 1.5 times the IQR of the hinges, points represent outliers beyond this range. When re-classifying recipients from 2006-2015 based on new Status 1-6 system, 6,255 (56%) Status 1A recipients met Status 3 critieria, 1,254 Status 2 critieria (11%) and 255 met Status 1 criteria (2%). A total of 3,462 (31%) of Status 1A recipients would have been downgraded to Status 4 because of violation of the cardiogenic shock requirement at the time of transplant. All 7,250 Status 1B recipients were re-assigned to Status 4. Two-hundred and twenty-eight (22%) low priority Status 2 (three-status system) would have been assigned the higher Status 4 because of restrictive cardiomyopathy, congenital heart disease, hypertrophic cardiomyopathy, or amyloidosis. The recipients are the same in both the three-status and six-status models, however, the increased number of tiers in the six-tier model allows for wider variation in estimated survival benefits. Compared to the three-status system, the within-status standard deviation of five-year survival benefit decreased from 5.5% to 4.9%. The majority (66%) of within-status variation in survival benefit was still attributable to centers, with the remaining 34% attributable to variation in transplant year, waiting time, donor quality and ischemic time.
#for eTable4
cand_thor <- haven::read_sas("SAF 2018 Q3/cand_thor.sas7bdat", NULL) %>%
haven::zap_formats() %>% haven::zap_labels()
init_lists <- to_model %>% group_by(PX_ID) %>%
filter(row_number()==1) %>%
ungroup() %>%
mutate(PX_ID = as.numeric(PX_ID)) %>%left_join(cand_thor, by = "PX_ID")
mean_age <- mean(init_lists$CAN_AGE_AT_LISTING)
pct_female <- 100*(init_lists %>% filter(CAN_GENDER == "F") %>% nrow())/(init_lists %>% nrow())
Of 29,172 candidates (mean age, 52; 26% women) listed at 113 centers.
#population estimates
##HR
pop_mean_hr <- comma_2(mean(pop_tx_hz$hr))
up_IQR_hr <- comma_2(quantile(pop_tx_hz$hr, probs = 0.75))
low_IQR_hr <- comma_2(quantile(pop_tx_hz$hr, probs = 0.25))
##five-year surival benefit
mean_benefit <- comma_2(100*mean(pop_tx_hz$benefit, na.rm = TRUE))
benefit_25 <- comma_2(100*quantile(pop_tx_hz$benefit, na.rm = TRUE, probs = 0.25))
benefit_75 <- comma_2(100*quantile(pop_tx_hz$benefit, na.rm = TRUE, probs = 0.75))
##five-year waitlist
mean_wait <- comma_2(100*mean(pop_tx_hz$marginal_waits, na.rm = TRUE))
wait_25 <- comma_2(100*quantile(pop_tx_hz$marginal_waits, na.rm = TRUE, probs = 0.25))
wait_75 <- comma_2(100*quantile(pop_tx_hz$marginal_waits, na.rm = TRUE, probs = 0.75))
##five-year post-transplant
mean_post <- comma_2(100*mean(pop_tx_hz$marginal_posts, na.rm = TRUE))
post_25 <- comma_2(100*quantile(pop_tx_hz$marginal_posts, na.rm = TRUE, probs = 0.25))
post_75 <- comma_2(100*quantile(pop_tx_hz$marginal_posts, na.rm = TRUE, probs = 0.75))
Heart transplantation was associated with a mean 76% decrease in mortality risk compared to waiting without transplant (hazard ratio 0.24, interquartile range [IQR] 0.14-0.35). This risk reduction generated a mean 44% (IQR 27-59%) absolute improvement in five-year survival, increasing from 33% (IQR 17- 51%) without transplant to 77% (IQR 74-80%) with transplant (full model results in Table S2).
#Calculate the within-status variation in the benefit under the three-status system
B_model <- lmer(100*benefit ~ stat_1b + stat_2 + (1|center), data = pop_tx_hz)
benefit_var <- as.data.frame(VarCorr(B_model))
center_var <- filter(benefit_var, grp == "center")$vcov
residual_var <- filter(benefit_var, grp == "Residual")$vcov
tot_var <- center_var + residual_var
sd_3 <- comma_2((tot_var)^(1/2))
pct_variation_center <- comma_2(100*center_var/tot_var)
pct_other <- comma_2(100*(1- center_var/tot_var))
#difference and 95% CI
diff_95ci <- function(x, y, w=NULL){
lm <- lm(y ~ x, weights = w)
diff <- comma_2(lm$coefficients[[2]])
confint(lm, level = 0.95)
low_ci <- comma_2(confint(lm, level = 0.95)[2,1])
up_ci <- comma_2(confint(lm, level = 0.95)[2,2])
if (lm$coefficients[[2]]>0){
return(paste0("+", diff, "; 95% CI, ", low_ci, " to ", up_ci))
}
paste0(diff, "; 95% CI, ", low_ci, " to ", up_ci)
}
diff_95ci_pct <- function(x, y, w=NULL){
lm <- lm(y ~ x, weights = w)
diff <- comma_2(100*lm$coefficients[[2]])
confint(lm, level = 0.95)
low_ci <- comma_2(100*confint(lm, level = 0.95)[2,1])
up_ci <- comma_2(100*confint(lm, level = 0.95)[2,2])
if (lm$coefficients[[2]]>0){
return(paste0("+", diff, "%; 95% CI, ", low_ci, "% to ", up_ci, "%"))
}
paste0(diff, "%; 95% CI, ", low_ci, "% to ", up_ci, "%")
}
high_low <- five_year_means %>% filter(sig_better ==1 | sig_worse ==1)
high_improvement <- comma(100*summary(lm(mean_five_year ~ sig_better,
data = high_low,
weights = high_low$num_tx))$coefficients[,1][[2]])
low_benefit <- five_year_means %>% filter(sig_worse ==1)
low_wait <- comma(100*weighted.mean(low_benefit$mean_five_year_wait,
low_benefit$num_tx))
low_post <- comma(100*weighted.mean(low_benefit$mean_five_year_post,
low_benefit$num_tx))
low_B <- comma(100*weighted.mean(low_benefit$mean_five_year,
low_benefit$num_tx))
high_benefit <- five_year_means %>% filter(sig_better ==1)
high_wait <- comma(100*weighted.mean(high_benefit$mean_five_year_wait,
high_benefit$num_tx))
high_post <- comma(100*weighted.mean(high_benefit$mean_five_year_post,
high_benefit$num_tx))
high_B <- comma(100*weighted.mean(high_benefit$mean_five_year,
high_benefit$num_tx))
post_high_low <- diff_95ci_pct(x = high_low$center_performance,
y = high_low$mean_five_year_post,
w = high_low$num_tx)
wait_high_low <- diff_95ci_pct(x = high_low$center_performance,
y = high_low$mean_five_year_wait,
w = high_low$num_tx)
surv_ben_high_low <- diff_95ci_pct(x = high_low$center_performance,
y = high_low$mean_five_year,
w = high_low$num_tx)
Adjusted for Status, candidate risk of death without transplant varied from 29% lower (relative HR 0.71) to 63% greater (relative HR 1.63) relative to the average center. The reduction in mortality risk from heart transplantation also varied significantly by center, from a maximum 29% greater reduction in mortality (relative HR 0.71) to a minimum of 68% lower reduction (relative HR 1.68) (eFigure 2). In combination, these center effects led to wide variation in the case-mix adjusted survival benefit of heart transplantation: the mean improvement in estimated five-year survival ranged from 30% to 55% by center, IQR [40% - 47%] (Figure 2). For every 10% decrease in expected candidate five-year waitlist survival without transplantation, there was a 6.2% increase in estimated survival benefit from heart transplant (95% CI 5.2% to 7.3%, Figure 3A). In contrast, there was no significant relationship between post-transplant survival and benefit (+1.5%, 95% CI -3.8% to 0.83%) (Figure 3B).
#Status at transplant by center performance
stat_by_benefit <- pop_tx_hz %>%
filter(tx ==1) %>%
left_join(five_year_means %>% select(center, center_performance)) %>%
filter(center_performance != "Average") %>%
mutate(center_performance = droplevels(center_performance),
six_status = factor(six_status,
labels = c("Status 1", "Status 2", "Status 3", "Status 4", "Status 6")))
pct_1a_High <- comma_2(100*mean(filter(stat_by_benefit, center_performance == "High")$stat_1a))
pct_1a_Low <- comma_2(100*mean(filter(stat_by_benefit, center_performance == "Low")$stat_1a))
p_1a_by_center <- comma_p(summary(glm(data =stat_by_benefit,
stat_1a ~ center_performance),
family = binomial, link = logit)$coefficients[2,4])
pct_1a_diff <- diff_95ci_pct(stat_by_benefit$center_performance, stat_by_benefit$stat_1a, w = NULL)
#Calculate overtreatment fraction amongst non-cardiogenic shock candidates
over_rates<- recips %>% left_join(five_year_means) %>%
filter(six_status > 3 & center_performance != "Average") %>%
mutate(over_1a = ifelse(three_status == "Status 1A", 1, 0))
over_low <- comma_2(100*mean(filter(over_rates,
center_performance == "Low" & six_status >3)$over_1a))
over_high <- comma_2(100*mean(filter(over_rates,
center_performance == "High" & six_status >3)$over_1a))
p_over <- comma_p(summary(glm(over_1a ~ center_performance,
data = over_rates,
family = binomial()))$coefficients[2,4])
over_diff <- diff_95ci_pct(x = over_rates$center_performance, y = over_rates$over_1a, w = NULL)
stat_by_benefit <- stat_by_benefit %>%
mutate(stat_just = fct_recode(stat_just,
"Status 1A (MCS complication)" = "Status 1A - (MCS complication)"))
#Status 1A characteristics
stat_1a_recip_w_tx_hr <- stat_by_benefit %>%
filter(three_status == "Status 1A") %>%
mutate(PX_ID = as.numeric(PX_ID)) %>% left_join(tx_hr, by = "PX_ID") %>%
mutate(Gender = factor(CAN_GENDER, levels = c("F", "M")),
Race = factor(CAN_RACE),
Race = fct_lump(Race, n = 3),
Race = fct_recode(Race,
"White" = "8",
"Black" = "16",
"Hispanic" = "2000",
"Other" = "Other"),
Diagnosis = case_when(
CAN_DGN>999 & CAN_DGN<1007 ~ "Dilated cardiomyopathy, non-ischemic",
CAN_DGN == 1007 | CAN_DGN ==1200 ~ "Ischemic cardiomyopathy",
CAN_DGN>1048 & CAN_DGN< 1100 ~ "Restrictive cardiomyopathy",
TRUE ~ "Other"
),
Diagnosis = factor(Diagnosis,
levels = c("Dilated cardiomyopathy, non-ischemic",
"Ischemic cardiomyopathy",
"Restrictive cardiomyopathy",
"Other")),
Diabetes = case_when(
CAN_DIAB_TY>1 & CAN_DIAB_TY<6 ~ "History of DM",
CAN_DIAB_TY ==1 ~ "Non-diabetic",
TRUE ~ "Unknown"
),
female_gfr = if_else(CAN_GENDER == "F", 0.742, 1),
black_gfr = if_else(Race == "Black", 1.21, 1),
eGFR = 175*((REC_CREAT)^(-1.154))*(REC_AGE_AT_TX^(-0.203))*female_gfr*black_gfr,
Renal_Function = case_when(
REC_DIAL == "Y" ~ "Dialysis",
eGFR >= 60 ~ "GFR >= 60 ml/min/1.73 m^2",
eGFR>= 30 ~ "GFR >= 30 & <60 ml/min/1.73 m^2",
eGFR < 30 ~ "GFR < 30 ml/min/1.73 m^2",
TRUE ~ "Unknown"
),
Renal_Function = if_else(is.na(Renal_Function)==TRUE, "Unknown", Renal_Function),
body_surface_area = 0.007184*(REC_HGT_CM)^(0.725)*REC_WGT_KG^(0.425),
Cardiac_Index = as.numeric(REC_CARDIAC_OUTPUT/body_surface_area),
Functional_Status = case_when(
REC_FUNCTN_STAT == 1 | (REC_FUNCTN_STAT>2069) ~"Limited Impairment, 10-30%",
(REC_FUNCTN_STAT>2039 & REC_FUNCTN_STAT<2061) ~ "Moderate Impairment, 40-60%",
(REC_FUNCTN_STAT>2000 & REC_FUNCTN_STAT<2031) ~ "Severe Impairment, 70-100%",
TRUE ~ "Unknown"
),
Functional_Status = ifelse(is.na(Functional_Status), "Unknown", Functional_Status),
severe_impairment = ifelse(Functional_Status == "Severe Impairment, 70-100%", 1, 0),
stat_1a_just = droplevels(stat_just),
acute_mcs = ifelse(stat_1a_just == "Status 1A (MCS for shock)", 1, 0),
lvad_comp = ifelse(stat_1a_just == "Status 1A (MCS complication)", 1, 0),
pcwp_15 = ifelse(REC_PCW_MEAN < 15, 1, 0),
pcwp_15 = ifelse(is.na(REC_PCW_MEAN), 0, pcwp_15),
blood_type = factor(
case_when(
CAN_ABO %in% c("A", "A1", "A2") ~ "A",
CAN_ABO %in% c("A1B", "A2B") ~ "AB",
TRUE ~ CAN_ABO)
),
payor = case_when(
REC_PRIMARY_PAY %in% c(3,4,13) ~ "Medicare",
REC_PRIMARY_PAY ==2 ~ "Medicaid",
REC_PRIMARY_PAY == 1 ~ "Private",
TRUE ~ "Other"
)
)
#PCWP < 15
pct_under_15_PCWP_high <- comma_2(100*mean(filter(stat_1a_recip_w_tx_hr,
center_performance == "High")$pcwp_15,
na.rm = TRUE))
pct_under_15_PCWP_low <- comma_2(100*mean(filter(stat_1a_recip_w_tx_hr,
center_performance == "Low")$pcwp_15,
na.rm = TRUE))
pcwp_15_diff <- diff_95ci_pct(x = stat_1a_recip_w_tx_hr$center_performance,
y = stat_1a_recip_w_tx_hr$pcwp_15)
#PCWP Mean
pct_REC_PCW_MEAN_High <- comma(mean(filter(stat_1a_recip_w_tx_hr,
center_performance == "High")$REC_PCW_MEAN,
na.rm = TRUE))
pct_REC_PCW_MEAN_Low <- comma(mean(filter(stat_1a_recip_w_tx_hr,
center_performance == "Low")$REC_PCW_MEAN,
na.rm = TRUE))
pcwp_mean_diff <- diff_95ci(x = stat_1a_recip_w_tx_hr$center_performance,
y = stat_1a_recip_w_tx_hr$REC_PCW_MEAN)
#functional impairment
pct_severe_impairment_High <- comma_2(100*mean(filter(stat_1a_recip_w_tx_hr,
center_performance == "High")$severe_impairment,
na.rm = TRUE))
pct_severe_impairment_Low <- comma_2(100*mean(filter(stat_1a_recip_w_tx_hr,
center_performance == "Low")$severe_impairment,
na.rm = TRUE))
severe_impair_diff <- diff_95ci_pct(x =stat_1a_recip_w_tx_hr$center_performance,
y =stat_1a_recip_w_tx_hr$severe_impairment)
#acute MCS use
pct_acute_mcs_High <- comma(100*mean(filter(stat_1a_recip_w_tx_hr,
center_performance == "High")$acute_mcs,
na.rm = TRUE))
pct_acute_mcs_Low <- comma(100*mean(filter(stat_1a_recip_w_tx_hr,
center_performance == "Low")$acute_mcs,
na.rm = TRUE))
acute_mcs_diff <- diff_95ci_pct(x = stat_1a_recip_w_tx_hr$center_performance,
y = stat_1a_recip_w_tx_hr$acute_mcs)
#LVAD complications
pct_lvad_comp_High <- comma_2(100*mean(filter(stat_1a_recip_w_tx_hr,
center_performance == "High")$lvad_comp,
na.rm = TRUE))
pct_lvad_comp_Low <- comma_2(100*mean(filter(stat_1a_recip_w_tx_hr,
center_performance == "Low")$lvad_comp,
na.rm = TRUE))
p_lvad_comp_by_center <- comma_p(summary(glm(data =stat_1a_recip_w_tx_hr,
lvad_comp ~ center_performance,
family = binomial))$coefficients[2,4])
lvad_comp_diff <- diff_95ci_pct(x = stat_1a_recip_w_tx_hr$center_performance,
y = stat_1a_recip_w_tx_hr$lvad_comp)
#center volume and benefit association
volume_benefit_diff <- diff_95ci(x = five_year_means$num_tx,
y = 100*100*five_year_means$mean_five_year)
volume_wait_diff <- diff_95ci(x = five_year_means$num_tx,
y = 100*100*five_year_means$mean_five_year_wait)
volume_post_diff <- diff_95ci(x = five_year_means$num_tx,
y = 100*100*five_year_means$mean_five_year_post)
lvad_rates <- to_model %>%
filter(tx ==1) %>%
left_join(center_re %>% select(center, center_performance)) %>%
filter(center_performance != "Average") %>%
mutate(high_ben = ifelse(center_performance == "High", 1, 0))
pct_lvad_high <- comma_2(100*mean(filter(lvad_rates,
center_performance == "High")$cf_lvad,
na.rm = TRUE))
pct_lvad_low <- comma_2(100*mean(filter(lvad_rates,
center_performance == "High")$cf_lvad,
na.rm = TRUE))
Compared to average, 31 centers (27%) had significantly higher survival benefit and 30 centers (27%) had significantly lower benefit. At five-years, Status-adjusted post-transplant survival was similar at high and low-benefit centers (77.6% vs. 77.1%, +0.5%; 95% CI, -1.3% to 2.3%). However, estimated candidate survival without transplant was lower at high-benefit centers (29% vs 39%, -10%; 95% CI, -12% to -8.1%), leading to a larger five-year survival benefit (48.6% vs 38%, +11%; 95% CI, 9.3% to 12%). Compared to low-benefit centers, high-benefit centers used Status 1A qualifying therapies less frequently (50% vs. 63%, -13%; 95% CI, -15% to -11%, Table S3) and were less likely to treat candidates who did not meet hemodynamic requirements for cardiogenic shock with Status 1A qualifying IABPs or high-dose inotropes (25% vs. 31%, -5.1%; 95% CI, -7.2% to -3%). At the time of transplant, Status 1A candidates at high benefit centers had higher pulmonary capillary wedge pressures (mean 20.1 vs. 18.9 mmHg, +1.2; 95% CI, 0.7 to 1.6 and percentage less than 15 mmHg 25% vs. 32%, -7.1%; 95% CI, -9.2% to -4.9%), and worse functional status (55% vs. 42% with severe impairment requiring continuous hospitalization, (+13%; 95% CI, 11% to 15%) than low-benefit centers. High-benefit centers transplanted more recipients for acute hemodynamic decompensation requiring mechanical support (37.6% vs. 30.5%, +7.2%; 95% CI, 4.9% to 9.4%) and used the “device-related complication” Status 1A justification significantly less often (20% vs. 37%, -17%; 95% CI, -19% to -15%) (Table S4). At high benefit centers, 35% of Status 1A recipients had a BTT LVAD at the time of transplant compared to35% of recipients at low benefit centers (-14%; 95% CI, -16% to -12%). Across all centers, transplant volume was weakly associated with post-transplant survival (+0.59; 95% CI, 0.074 to 1.1), however was not associated with either candidate urgency (+0.49; 95% CI, -0.29 to 1.3) or survival benefit (+0.099; 95% CI, -0.56 to 0.75) (Figure S4).
#status switch table
status_switch <- pop_tx_hz %>% left_join(pop_tx_hz_6_stat) %>%
group_by(six_status, three_status) %>% count() %>%
spread(key = three_status, value = n) %>%
mutate_all(funs(ifelse(is.na(.)==TRUE, 0, .)))
tot_1as <- sum(status_switch$"Status 1A")
tot_1bs <- sum(status_switch$"Status 1B")
tot_2s <- sum(status_switch$"Status 2")
Under the prior three-status heart allocation system, the most common Status at transplant was 1A (N =11,227, 57%) and heart transplantation for a Status 1A candidate conferred a mean five-year survival benefit associated with 58%, (IQR 54- 62%). Status 1B was the next most common (N = 7,250, 37%) with a mean five-year survival benefit associated with transplant of 27%, (IQR 23- 31%). There were only 1,338 candidates at Status 2 at the time of transplant (7%) with mean five-year survival benefit associated with transplant of 14%, (IQR 10- 17%) (Figure 4A). Overall, within a given priority Status at transplantation, the standard deviation in five-year survival benefit was 5.5% with the majority (76%) of the variation amongst recipients attributable to center-level effects and the remaining 24% attributable to variation in transplant year, waiting time, donor quality and ischemic time. When retrospectively re-classifying recipients from 2006-2015 based on new Status 1-6 system, 6,255 Status 1A recipients met Status 3 criteria (56%), 1,254 met Status 2 criteria (11%), and 255 met Status 1 criteria (2%). A total of 3,462 (31%) of Status 1A recipients would have been downgraded to Status 4 because of a violation of the cardiogenic shock requirement at the time of transplant. All 7,250 Status 1B recipients were re-assigned to Status 4. Two-hundred and twenty eight (22%) low priority Status 2 (three-status system) recipients would have been assigned to the higher Status 4 because of restrictive cardiomyopathy, congenital heart disease, hypertrophic cardiomyopathy, or amyloidosis. (eTable5). Estimated five-year survival benefit associated with transplant ranged from of 68% (Status 1) to 14% (Status 6) (Figure 4B). Compared to the three-status system, the within-status standard deviation of five-year survival benefit decreased from 5.5% to 4.9% but the majority (66%) was still attributable to centers. The between-center variance in the survival benefit of transplant on the log hazard ratio scale transplant decreased from 0.161 to 0.115, corresponding to a reduction in the between-center variance of absolute five-year survival benefit from 23% to 16%. (eTable6).
#Variation in log hazard of transplant in three and six-status models
var_tx_3_stat <- VarCorr(three_status)$center[2,2]
var_tx_6_stat <- VarCorr(six_status)$center[2,2]
#within-status variation in the five-year survival benefit under the six-status system
B_model_6 <- lmer(100*benefit_6 ~ factor(six_status) + (1|center), data = pop_tx_hz_6_stat)
benefit_var_6 <- as.data.frame(VarCorr(B_model_6))
center_var_6 <- filter(benefit_var_6, grp == "center")$vcov
residual_var_6 <- filter(benefit_var_6, grp == "Residual")$vcov
tot_var_6 <- center_var_6 + residual_var_6
sd_6 <- comma_2((tot_var_6)^(1/2))
pct_variation_center_6 <- comma_2(100*center_var_6/tot_var_6)
pct_other_6 <- comma_2(100*(1- center_var_6/tot_var_6))
When re-classifying recipients from 2006-2015 based on new Status 1-6 system, the majority (56%) of Status 1A recipients would have been Status 3, with a smaller group qualifying for Status 2 (11%) and very few meeting Status 1 criteria (2%). A substantial portion (31%) of Status 1A candidates would have been downgraded to Status 4 because of violation of the cardiogenic shock requirement. A small number of recipients (N = 288) would have been assigned a higher status because of restrictive cardiomyopathy (Table S5 for details). Five-year survival benefit from transplant ranged from of 68% (Status 1) to 14% (Status 6) (Figure 5B).
Compared to the three-status system, the within-status standard deviation of five-year survival benefit decreased from 5.5% to 4.9% but the majority (66%) was still attributable to centers. The between-center variance in the survival benefit of transplant on the log hazard scale decreased from 0.161 to 0.115 on log hazard ratio scale. The between-center variance in absolute 5-year survival benefit decreased from 23% to 16%. (Table S6).
Current and future adult heart allocation. Schematic depiction of the shift from the current adult heart allocation system to the modified system. * Cardiogenic shock requirement applies. (Constructed with permission directly from the policy details in Organ Procurement and Transplantation Network1). Complete detail about each of the Status criteria are available at https://optn.transplant.hrsa.gov/media/1200/optn_policies.pdf#nameddest=Policy_06. Dischargeable vs. non-dischargeable ventricular support devices were distinguished via a list provided by the SRTR.
#Survival benefit from transplantation for a Status 1A patient
#transplanted at a center with a high benefit from transplant
top_c <- center_re %>% arrange(-benefit) %>% filter(row_number()==1)
top_c_tx_hr <- top_c$tx_exp
top_c_wt_risk <- top_c$intercept_exp
top_c <- top_c$center
best_ctr_recips <- unique(to_model %>% filter(center == top_c) %>%
select(PX_ID)) %>%
left_join(to_model %>% select(PX_ID, tx, stat_1a, t_1, era_tx) %>%
group_by(PX_ID, tx) %>%
filter(tx ==1 & row_number() ==1)) %>%
filter(is.na(tx) == FALSE)
ex_pt_top <- best_ctr_recips[[156, "PX_ID"]]
t_tx <- best_ctr_recips[[156, "t_1"]]
hz_good <- exp(hr_tx_ij(model = three_status,
df = to_model, i = top_c, j = ex_pt_top))
max_fup <- max(basehz$time)/365
ex_benefit <- ben_wait_tx(df = to_model,
base_hz = basehz,
model = three_status, j = ex_pt_top, i = top_c)
fiveyr_good <- comma_2(100*(ex_benefit[[2]]-ex_benefit[[1]]))
fiveyr_good_wait <- comma_2(100*ex_benefit[[1]])
fiveyr_good_tx <- comma_2(100*ex_benefit[[2]])
good_plot <- b_ij(df = to_model, base_hz = basehz, model = three_status, j = ex_pt_top, i = top_c) +
scale_x_continuous(limits = c(0, 5*365), breaks= seq(0, t_tx + 5*365, 365)) +
scale_y_continuous(breaks= seq(0, 1, 0.2)) +
theme_economist_white() +
theme(axis.title = element_text(),
axis.title.y = element_text(),
legend.position="bottom",
legend.key = element_rect(fill = "white"),
legend.background = element_rect(fill = "white"),
legend.text = element_text(size = 12))
### Figure 2B: Survival benefit from transplantation for a "Status 1A" patient
#transplanted at a center with a lower benefit from transplant
worst_c <- center_re %>% arrange(benefit) %>% filter(row_number()==1)
worst_c_tx_hr <- worst_c$tx_exp
worst_c_wt_risk <- worst_c$intercept_exp
worst_c <- worst_c$center
worst_c_recips <- unique(to_model %>%
filter(center == worst_c) %>%
select(PX_ID)) %>%
left_join(to_model %>% select(PX_ID, tx, stat_1a, t_1, era_tx) %>%
group_by(PX_ID, tx) %>%
filter(tx ==1 & row_number() ==1)) %>%
filter(is.na(tx) == FALSE)
ex_worst <- worst_c_recips[[30,"PX_ID"]]
t_tx_worst <- worst_c_recips[[30,"t_1"]]
ex_benefit_bad <- ben_wait_tx(df = to_model,
base_hz = basehz,
model = three_status,
i = worst_c,
j = ex_worst)
fiveyr_bad <- comma_2(100*(ex_benefit_bad[[2]]-ex_benefit_bad[[1]]))
fiveyr_bad_wait <- comma_2(100*ex_benefit_bad[[1]])
fiveyr_bad_tx <- comma_2(100*ex_benefit_bad[[2]])
hz_bad <- exp(hr_tx_ij(model = three_status, df = to_model, i = worst_c, j = ex_worst))
bad_plot <- b_ij(df = to_model,
base_hz = basehz,
model = three_status,
i = worst_c, j = ex_worst)+
scale_x_continuous(limits = c(0, 5*365),
breaks= seq(0, t_tx_worst + 5*365, 365)) +
scale_y_continuous(breaks= seq(0, 1, 0.2)) +
theme_economist_white() +
theme(axis.title = element_text(),
axis.title.y = element_text(), legend.position="right")
#combine plots
g_legend<-function(a.gplot){
tmp <- ggplot_gtable(ggplot_build(a.gplot))
leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box")
legend <- tmp$grobs[[leg]]
return(legend)}
mylegend <- g_legend(good_plot)
fig_2 <- grid.arrange(
arrangeGrob(good_plot + theme(legend.position="none"),
bad_plot + theme(legend.position="none",
axis.title.y = element_blank(),
axis.text.y = element_blank()),
nrow=1, widths = c(1, 0.9)),
mylegend,
nrow=2,heights=c(10, 2.2))
ggsave("efig_2.pdf", plot = fig_2)
Two example survival benefit predictions incorporating the timing of transplantation, candidate status, donor quality, ischemic time, and year of transplant.
The high-benefit center (Panel A) prioritizes Status 1A candidates at higher risk of wait-list mortality and transplantation at this center lowers risk more than average. Specifically, candidates listed by this center have a risk of death 63% higher than the average center (relative HR 1.63) and transplant lowers the risk of mortality 29% more than at the average center (HR 0.71). The combination of these two effects leads to a large 5-year survival benefit, improving survival (conditional upon reaching transplant) from 6.1% to 75% (+69% absolute gain) for this particular Status 1A patient transplanted after waiting 135 days.
In contrast, the patients transplanted at the low-benefit center (Panel B) are at lower risk of dying on the waitlist at baseline and transplantation at this center lowers risk less than average. Specifically, candidates listed by this center have a risk of death 26% lower than the average center (HR 0.74) and transplant lowers the risk of mortality 63% less than at the average center (HR 1.63). Despite the fact that the candidate has the same priority Status as the candidate in Panel A, the combination of these two effects leads to a smaller 5-year survival benefit, improving survival from 28% to 75% (+46% absolute gain) for transplanted after waiting 162 days.
coxme resultsprint(three_status, digits =2)
Cox mixed-effects model fit by maximum likelihood
Data: to_model
events, n = 11058, 163240
Iterations= 30 184
NULL Integrated Fitted
Log-likelihood -107600 -105890.7 -105702
Chisq df p AIC BIC
Integrated loglik 3419 10 0 3399 3326
Penalized loglik 3796 123 0 3550 2653
Model: f_coxme
Fixed coefficients
coef exp(coef) se(coef) z p
tx -1.953 0.14 0.045 -43.7 0.000
stat_2 -1.355 0.26 0.036 -37.1 0.000
stat_1b -0.942 0.39 0.032 -29.8 0.000
tx_risky_don 0.260 1.30 0.028 9.3 0.000
era_tx -0.054 0.95 0.030 -1.8 0.077
tx:stat_2 1.321 3.75 0.064 20.7 0.000
tx:stat_1b 0.904 2.47 0.044 20.7 0.000
Random effects
Group Variable Std Dev Variance Corr
center Intercept 0.224 0.050 -0.509
tx 0.238 0.057
table_fe_results <- function(me_cox_1){
fe_results_table <- tibble(covariate = character(),
`log hazard ratio` = character(),
`95% CI` = character())
fe <- fixef(me_cox_1)
vc <- vcov(me_cox_1)
for (i in seq_along(fe)){
name <- names(fe[i])
hr <- comma(unname(fe[i]))
se <- sqrt(vc[i,i])
ci <- paste0("(", comma(unname(fe[i]) - 1.96*se),
", ",
comma(unname(fe[i]) + 1.96*se), ")")
fe_results_table <- rbind(fe_results_table,
tibble(covariate =name,
`log hazard ratio` = hr,
`95% CI` = ci))
}
return(fe_results_table)
}
table_fe_results(three_status)
#write_csv(table_fe_results(three_status), "three_status_fixed.csv")
to_plot <- plot_tx_effects(three_status) %>%
mutate(`Patient Status` = factor(covariate,
levels = c("tx",
"tx:stat_1b",
"tx:stat_2",
"tx_risky_don",
"era_tx"),
labels = c("Status 1A",
"Status 1B",
"Status 2",
"High risk donor (Status 1A recipient)",
"Before 2010 (Status 1A recipient)"))) %>%
arrange(factor(covariate,
levels = c("tx","tx:stat_1b", "tx:stat_2", "tx_risky_don","era_tx")))
ggplot(data = to_plot %>%
filter(!covariate %in% c("tx_risky_don","era_tx")),
aes(y = `Patient Status`, x = haz)) +
geom_point() +
geom_errorbarh(aes(xmin = low_ci, xmax = up_ci)) +
theme_fivethirtyeight() +
theme(axis.title = element_text(), axis.title.y = element_blank()) +
scale_x_continuous(limits = c(0,1),
breaks = c(0, 0.25, 0.5, 0.75, 1),
name = "Benefit heart transplantation (HR for Death)" )
tx_est <- comma_2(filter(to_plot, covariate == "tx")$haz)
tx_low <-comma_2(filter(to_plot, covariate == "tx")$low_ci)
tx_up <- comma_2(filter(to_plot, covariate == "tx")$up_ci)
stat_1b_est <- comma_2(filter(to_plot, covariate == "tx:stat_1b")$haz)
stat_1b_low <-comma_2(filter(to_plot, covariate == "tx:stat_1b")$low_ci)
stat_1b_up <- comma_2(filter(to_plot, covariate == "tx:stat_1b")$up_ci)
stat_2_est <- comma_2(filter(to_plot, covariate == "tx:stat_2")$haz)
stat_2_low <-comma_2(filter(to_plot, covariate == "tx:stat_2")$low_ci)
stat_2_up <- comma_2(filter(to_plot, covariate == "tx:stat_2")$up_ci)
ggsave("HR_transplant_model_1.pdf")
Hazard ratios of transplantation in the old 3-status system at transplantation as estimated by mixed-effects model, adjusted for donor risk, ischemic time, and year of transplant (see supplement for full model results).
Status 1A candidates had a hazard ratio of transplant 0.14 (95% CI 0.13-0.15), Status 1B candidates had a HR of 0.35 (95% CI 0.32-0.38) and Status 2 candidates had a HR of 0.53 (95% CI 0.47-0.61).
ggplot(data = basehz, aes(x=time/365, y = hz_t)) +
geom_smooth() +
labs(y = "Hazard") +
scale_x_continuous(name = "Years after listing", limits = c(0, 12), breaks = seq(0, 12, 2))
ggsave("base_hazard.pdf")
This baseline hazard function used to calculating the 5-year benefit was estimated non-parametrically using a kaplan-meier estimate. Immediately after listing for heart transplantation, candidates have increased hazard of death. More than 5 years after listing the hazard begins to rise again due to aging.
##Figure S3: Distribution of Survival Benefit at the Best and Worst Center, 2006-2015
best_center <- paste0(as.character(
five_year_means[five_year_means$mean_five_year ==max(five_year_means$mean_five_year),]$center), ".x")
worst_center <- paste0(as.character(
five_year_means[five_year_means$mean_five_year ==min(five_year_means$mean_five_year),]$center), ".x")
to_plot <- ben_five_year_only %>% select(best_center, worst_center) %>%
gather(key = center, value = five_year_benefit) %>%
mutate(ctr = if_else(center == best_center, "Best Center", "Worst Center"))
best_ctr_mean <- comma_2(
100*mean(filter(to_plot, ctr == "Best Center")$five_year_benefit, na.rm = TRUE))
worst_ctr_mean <- comma_2(
100*mean(filter(to_plot, ctr == "Worst Center")$five_year_benefit, na.rm =TRUE))
ggplot(to_plot, aes(x = five_year_benefit*100, fill = ctr)) +
geom_density(alpha = 0.5) +
scale_x_continuous(limits = c(-25, 100), name = "5-year survival benefit (%)") +
scale_fill_manual(values = c("blue", "red")) +
guides(fill = guide_legend(""))
ggsave("fig_S4.pdf")
Distribution of 5-year survival benefit if the best and worst center hypothetically had transplanted all recipients during the study time period. The best center had a mean five year survival benefit of 55% compared to only 30% at the lowest benefit center.
tangram::rmd(tangram::tangram("center_performance ~ three_status + six_status", stat_by_benefit, digits = 3))
| N | Low | High | Test Statistic | |
|---|---|---|---|---|
| (N=5430) | (N=6878) | |||
| three_status | 12308 | X22=238.05, P=<0.001 | ||
| Status 1A | 0.632 3430/5430 | 0.501 3447/6878 | ||
| Status 1B | 0.331 1800/5430 | 0.423 2908/6878 | ||
| Status 2 | 0.037 200/5430 | 0.076 523/6878 | ||
| six_status | 12308 | X24=285.82, P=<0.001 | ||
| Status 1 | 0.012 66/5430 | 0.016 110/6878 | ||
| Status 2 | 0.075 407/5430 | 0.051 349/6878 | ||
| Status 3 | 0.383 2079/5430 | 0.265 1820/6878 | ||
| Status 4 | 0.502 2725/5430 | 0.613 4214/6878 | ||
| Status 6 | 0.028 153/5430 | 0.056 385/6878 |
#Label Variables when necessary
Hmisc::label(stat_1a_recip_w_tx_hr$Cardiac_Index) <- "Cardiac Index (L/min/m^2)"
Hmisc::label(stat_1a_recip_w_tx_hr$REC_PCW_MEAN) <- "Mean Pulmonary Capillary Wedge Pressure (mmHg)"
Hmisc::label(stat_1a_recip_w_tx_hr$stat_1a_just) <- "Status 1A Justification at Transplant"
tangram::rmd(tangram::tangram("center_performance ~ benefit + marginal_posts + marginal_waits + t_1 + REC_AGE_AT_TX + REC_BMI + Gender + Race + Diagnosis + Diabetes + Renal_Function + Functional_Status + Cardiac_Index + REC_PCW_MEAN + stat_1a_just + payor", data = stat_1a_recip_w_tx_hr))
| N | Low | High | Test Statistic | |
|---|---|---|---|---|
| (N=3430) | (N=3447) | |||
| benefit | 6836 | 0.505 0.534 0.558 | 0.606 0.635 0.651 | F1,6834=17466.59, P=<0.001 |
| marginal_posts | 6836 | 0.740 0.774 0.806 | 0.740 0.776 0.799 | F1,6834=0.40, P=0.528 |
| marginal_waits | 6836 | 0.220 0.241 0.262 | 0.120 0.138 0.167 | F1,6834=13817.14, P=<0.001 |
| t_1 | 6877 | 43 138 374 | 16 56 213 | F1,6875=337.56, P=<0.001 |
| Age at TX | 6877 | 45 55 62 | 45 55 62 | F1,6875=0.61, P=0.436 |
| BMI: | 6866 | 23.8 27.2 31.0 | 23.3 26.5 30.1 | F1,6864=30.68, P=<0.001 |
| Gender : M | 6877 | 0.788 2703/3430 | 0.774 2668/3447 | X21=1.98, P=0.159 |
| Race | 6877 | X23=28.01, P=<0.001 | ||
| White | 0.655 2246/3430 | 0.634 2186/3447 | ||
| Black | 0.226 775/3430 | 0.206 709/3447 | ||
| Hispanic | 0.084 288/3430 | 0.104 359/3447 | ||
| Other | 0.035 121/3430 | 0.056 193/3447 | ||
| Diagnosis | 6877 | X23=7.13, P=0.068 | ||
| Dilated cardiomyopathy, non-ischemic | 0.456 1563/3430 | 0.457 1577/3447 | ||
| Ischemic cardiomyopathy | 0.359 1232/3430 | 0.346 1191/3447 | ||
| Restrictive cardiomyopathy | 0.107 368/3430 | 0.102 352/3447 | ||
| Other | 0.078 267/3430 | 0.095 327/3447 | ||
| Diabetes | 6877 | X22=2.66, P=0.264 | ||
| History of DM | 0.291 997/3430 | 0.274 945/3447 | ||
| Non-diabetic | 0.708 2429/3430 | 0.724 2496/3447 | ||
| Unknown | 0.001 4/3430 | 0.002 6/3447 | ||
| Renal_Function | 6877 | X24=36.83, P=<0.001 | ||
| Dialysis | 0.032 110/3430 | 0.031 107/3447 | ||
| GFR < 30 ml/min/1.73 m^2 | 0.025 87/3430 | 0.038 131/3447 | ||
| GFR >= 30 & <60 ml/min/1.73 m^2 | 0.361 1238/3430 | 0.386 1329/3447 | ||
| GFR >= 60 ml/min/1.73 m^2 | 0.582 1995/3430 | 0.540 1860/3447 | ||
| Unknown | 0.000 0/3430 | 0.006 20/3447 | ||
| Functional_Status | 6877 | X23=134.07, P=<0.001 | ||
| Limited Impairment, 10-30% | 0.264 907/3430 | 0.172 592/3447 | ||
| Moderate Impairment, 40-60% | 0.268 920/3430 | 0.237 817/3447 | ||
| Severe Impairment, 70-100% | 0.424 1455/3430 | 0.553 1907/3447 | ||
| Unknown | 0.043 148/3430 | 0.038 131/3447 | ||
| Cardiac Index | 6458 | 1.80 2.22 2.69 | 1.76 2.16 2.64 | F1,6456=9.25, P=0.002 |
| Mean Pulmonary Capillary Wedge Pressure | 6181 | 12.0 18.0 25.0 | 13.0 20.0 26.0 | F1,6179=30.51, P=<0.001 |
| Status 1A Justification at Transplant | 6877 | X25=247.04, P=<0.001 | ||
| Status 1A - (Exception) | 0.072 246/3430 | 0.087 299/3447 | ||
| Status 1A - (High dose inotropes) | 0.250 858/3430 | 0.333 1149/3447 | ||
| Status 1A (MCS complication) | 0.369 1266/3430 | 0.200 690/3447 | ||
| Status 1A - (Mechanical ventilation) | 0.003 12/3430 | 0.003 12/3447 | ||
| Status 1A (MCS for shock) | 0.305 1045/3430 | 0.376 1297/3447 | ||
| Status 1A (no justification listed) | 0.001 3/3430 | 0.000 0/3447 | ||
| payor | 6877 | X23=76.76, P=<0.001 | ||
| Medicaid | 0.124 426/3430 | 0.128 442/3447 | ||
| Medicare | 0.365 1251/3430 | 0.301 1039/3447 | ||
| Other | 0.055 188/3430 | 0.028 96/3447 | ||
| Private | 0.456 1565/3430 | 0.543 1870/3447 |
For continuous variables, median and IQRs are presented and group comparison testing p-values obtained using Kruskal-Wallis. For categorical variables, number (%) are presented and group comparison testing p-values calculated with Pearson Chi-squared test. GFR = Glomerular Filtration Rate, calculated by the Modification in Diet in Renal Disease (MDRD) formula. Body surface area for cardiac index calculated with method of Du Bois. Functional status is on the Karnofsky scale.
tangram::rmd(tangram::tangram(data = stat_by_benefit, three_status~ six_status))
| N | Status 1A | Status 1B | Status 2 | Test Statistic | |
|---|---|---|---|---|---|
| (N=6877) | (N=4708) | (N=723) | |||
| six_status | 12308 | X28=15011.30, P=<0.001 | |||
| Status 1 | 0.026 176/6877 | 0.000 0/4708 | 0.000 0/723 | ||
| Status 2 | 0.110 756/6877 | 0.000 0/4708 | 0.000 0/723 | ||
| Status 3 | 0.567 3899/6877 | 0.000 0/4708 | 0.000 0/723 | ||
| Status 4 | 0.298 2046/6877 | 1.000 4708/4708 | 0.256 185/723 | ||
| Status 6 | 0.000 0/6877 | 0.000 0/4708 | 0.744 538/723 |
Classification of recipients from 2006-2015 based on new Status 1-6 system. The majority ( 56%) of Status 1A recipients were coded as Status 3. A substantial portion (31%) of Status 1A candidates would have been downgraded to Status 4 because of violation of the cardiogenic shock requirement. A small number of recipients 288 were assigned a higher status because of restrictive cardiomyopathy
##Figure S4: Association of center transplant volume with wait-list survival, post-transplant survival, and survival benefit from heart transplantation.
to_plot <- five_year_means %>%
select(num_tx, mean_five_year, mean_five_year_post, mean_five_year_wait) %>%
gather(key, value, -num_tx) %>%
mutate( key = case_when(
key == "mean_five_year" ~ "survival benefit",
key == "mean_five_year_wait" ~ "waitlist survival",
key == "mean_five_year_post" ~ "post-transplant survival"
)
)
ggplot(to_plot, aes(x = num_tx, y = 100*value, color = key, group = key)) +
geom_point( size = 3) + geom_smooth( method = "lm") +
scale_color_brewer(palette = "Set1", direction = -1 ) +
labs(y = "Five-year estimate (%)",
x = "Number of Transplants performed during study time period",
color = "") +
ylim(0, 100)
ggsave("Fig_S4.pdf")
Across all centers, transplant volume was weakly associated with post-transplant survival (+0.59; 95% CI, 0.074 to 1.1), however was not associated with either candidate urgency (+0.49; 95% CI, -0.29 to 1.3) or survival benefit (+0.099; 95% CI, -0.56 to 0.75) (Figure S4).
coxme resultsprint(six_status, digits = 2)
Cox mixed-effects model fit by maximum likelihood
Data: to_model
events, n = 11058, 163240
Iterations= 23 211
NULL Integrated Fitted
Log-likelihood -107600 -105838.3 -105664.7
Chisq df p AIC BIC
Integrated loglik 3523 14 0 3495 3393
Penalized loglik 3871 125 0 3621 2707
Model: f_coxme
Fixed coefficients
coef exp(coef) se(coef) z p
tx -2.924 0.054 0.139 -21.0 0.0e+00
status_2 -0.869 0.419 0.089 -9.7 0.0e+00
status_3 -1.526 0.217 0.083 -18.3 0.0e+00
status_4 -2.169 0.114 0.079 -27.6 0.0e+00
status_6 -2.614 0.073 0.082 -31.8 0.0e+00
tx_risky_don 0.262 1.300 0.028 9.4 0.0e+00
era_tx -0.053 0.949 0.030 -1.7 8.3e-02
tx:status_2 0.737 2.090 0.152 4.8 1.3e-06
tx:status_3 1.250 3.489 0.142 8.8 0.0e+00
tx:status_4 1.851 6.368 0.138 13.4 0.0e+00
tx:status_6 2.287 9.841 0.149 15.3 0.0e+00
Random effects
Group Variable Std Dev Variance Corr
center Intercept 0.195 0.038 -0.386
tx 0.212 0.045
#write_csv(table_fe_results(six_status), "six_status_fixed.csv")
table_fe_results(six_status)
to_plot <- plot_tx_effects(six_status) %>%
mutate(`Patient Status` = factor(covariate,
levels = c("tx",
"tx:status_2",
"tx:status_3",
"tx:status_4",
"tx:status_6",
"tx_risky_don","era_tx"),
labels = c("Status 1",
"Status 2",
"Status 3",
"Status 4",
"Status 6",
"High risk donor (Status 1 recipient)",
"Before 2010 (Status 1 recipient)")))
ggplot(data = to_plot %>%
filter(!(covariate %in% c("tx_risky_don","era_tx"))),
aes(y = `Patient Status`, x = haz)) +
geom_point() +
geom_errorbarh(aes(xmin = low_ci, xmax = up_ci)) +
scale_x_continuous(limits = c(0,0.75), breaks = c(0, 0.25, 0.5, 0.75),
name = "Hazard Ratio for death relative to waiting") +
theme_fivethirtyeight() +
theme(axis.title = element_text(),
axis.title.y = element_blank(),
axis.title.x = element_text(vjust = 0.1))
stat_list <- c("status_1", "status_2", "status_3", "status_4", "status_6")
tx_estimates <- tibble(status = stat_list) %>%
mutate(
hazard = "",
low_ci = "",
up_ci = ""
)
i <- 1
for (s in stat_list) {
if (s == "status_1"){
int <- "tx"
}
else{
int <- paste0("tx:", s)
}
tx_estimates$hazard[i] <-comma_2(filter(to_plot, covariate == int)$haz)
tx_estimates$low_ci[i] <-comma_2(filter(to_plot, covariate == int)$low_ci)
tx_estimates$up_ci[i] <-comma_2(filter(to_plot, covariate == int)$up_ci)
i <- i + 1
}
ggsave("HR_transplant_six_status.pdf")
Hazard ratios of transplantation by 6-Status Justification at transplantation as estimated by mixed-effects model, adjusted for donor risk, ischemic time, and year of transplant (see supplement for full model results).
Hazard ratios of transplantation conditional Status at transplantation as estimated by mixed-effects model (see supplement for full model results). Status 1 candidates had hazard ratio of transplant of 0.054 (95% CI 0.041-0.071). Compared to this group, Status 2-6 candidates had significantly lower benefit from heart transplanation.
Status 2 had a hazard ratio of 0.11 (95% CI 0.096-0.13), Status 3 0.19 (95% CI 0.17-0.21), Status 4 0.34 (95% CI 0.32-0.37), and Status 6 0.53 (95% CI 0.46-0.61)
Here we fit the model
\(h(t) = h_0(t)exp(center + status + tx*(center + status + donrisk + txYear)\)
num_tx_by_center <- to_model %>%
filter(dead ==1) %>%
group_by(center) %>%
count()
#center 602 is close to the "mean" center
for_fe <- to_model %>%
left_join(num_tx_by_center, by = "center") %>%
filter(n>20) %>%
mutate(three_status = factor(three_status),
center = factor(center),
center = relevel(center, ref = "602"))
full_fe <- "Surv(t_1, t_2, dead) ~ tx*(three_status + center) + tx_risky_don + era_tx"
full_fe_model <- coxph(as.formula(full_fe), for_fe)
centers_exlcuded <- length(unique(to_model$center)) - length(levels(for_fe$center))
This model treats center as a fixed effect. To get this model to converge, we excluded 9 centers with less than 20 deaths (before or after transplant)
print(full_fe_model, digits = 2)
Call:
coxph(formula = as.formula(full_fe), data = for_fe)
coef exp(coef) se(coef) z p
tx -1.7846 0.1679 0.2160 -8.3 <2e-16
three_statusStatus 1B -0.9490 0.3871 0.0320 -29.6 <2e-16
three_statusStatus 2 -1.3824 0.2510 0.0371 -37.3 <2e-16
center7 -0.0589 0.9428 0.1853 -0.3 0.751
center8 -0.4979 0.6078 0.2253 -2.2 0.027
center18 -0.4205 0.6567 0.2099 -2.0 0.045
center25 -0.1087 0.8970 0.2166 -0.5 0.616
center33 -0.2485 0.7800 0.1753 -1.4 0.156
center47 0.1439 1.1548 0.2626 0.5 0.584
center55 0.4539 1.5745 0.2510 1.8 0.071
center63 -0.2472 0.7810 0.2720 -0.9 0.363
center64 0.1690 1.1841 0.2253 0.8 0.453
center65 -0.4987 0.6073 0.2625 -1.9 0.057
center67 -0.1665 0.8466 0.2125 -0.8 0.433
center72 0.0993 1.1044 0.1916 0.5 0.604
center78 -0.0439 0.9570 0.1759 -0.2 0.803
center79 -0.1057 0.8997 0.2104 -0.5 0.616
center91 -0.1695 0.8441 0.2152 -0.8 0.431
center94 -0.6451 0.5246 0.2511 -2.6 0.010
center96 -0.2542 0.7755 0.2185 -1.2 0.245
center108 -0.2163 0.8055 0.1856 -1.2 0.244
center119 -0.0430 0.9580 0.1894 -0.2 0.821
center128 -0.5579 0.5724 0.2086 -2.7 0.007
center131 -0.4926 0.6110 0.1951 -2.5 0.012
center136 -0.3042 0.7377 0.1843 -1.7 0.099
center141 -0.0155 0.9846 0.1747 -0.1 0.929
center147 -0.1943 0.8234 0.2512 -0.8 0.439
center149 -0.2963 0.7435 0.2719 -1.1 0.276
center165 -0.4731 0.6230 0.2389 -2.0 0.048
center171 -0.7144 0.4895 0.2201 -3.2 0.001
center178 -0.5852 0.5570 0.2024 -2.9 0.004
center183 -0.3760 0.6866 0.2066 -1.8 0.069
center185 0.2324 1.2616 0.2446 1.0 0.342
center190 0.0522 1.0535 0.1979 0.3 0.792
center195 -0.1341 0.8745 0.2339 -0.6 0.566
center199 -0.1887 0.8281 0.2719 -0.7 0.488
center201 0.0038 1.0038 0.1958 0.0 0.984
center210 -0.0412 0.9596 0.1832 -0.2 0.822
center213 -0.5390 0.5833 0.2672 -2.0 0.044
center221 0.0137 1.0138 0.1739 0.1 0.937
center234 0.4450 1.5605 0.3046 1.5 0.144
center298 0.0989 1.1040 0.1815 0.5 0.586
center301 -0.5336 0.5865 0.1779 -3.0 0.003
center303 -0.4471 0.6395 0.1910 -2.3 0.019
center307 -0.3072 0.7355 0.1979 -1.6 0.121
center313 0.3507 1.4201 0.1843 1.9 0.057
center324 -0.3478 0.7062 0.2170 -1.6 0.109
center337 -0.2016 0.8175 0.1760 -1.1 0.252
center338 -0.2270 0.7969 0.2416 -0.9 0.347
center348 -0.1250 0.8825 0.1739 -0.7 0.472
center350 -0.4058 0.6664 0.1812 -2.2 0.025
center351 0.0669 1.0692 0.1679 0.4 0.690
center359 -0.7082 0.4925 0.2445 -2.9 0.004
center372 0.2889 1.3349 0.2183 1.3 0.186
center377 -0.6394 0.5276 0.3134 -2.0 0.041
center380 -0.2149 0.8066 0.2219 -1.0 0.333
center382 -0.0611 0.9407 0.1735 -0.4 0.725
center388 -0.0876 0.9162 0.2124 -0.4 0.680
center402 -0.7369 0.4786 0.2477 -3.0 0.003
center404 -0.1221 0.8851 0.1953 -0.6 0.532
center408 -0.4545 0.6348 0.2586 -1.8 0.079
center434 -0.3670 0.6928 0.1617 -2.3 0.023
center442 -0.2675 0.7653 0.1798 -1.5 0.137
center445 -0.3730 0.6887 0.2168 -1.7 0.085
center446 -0.4926 0.6110 0.1966 -2.5 0.012
center459 -0.1680 0.8453 0.2256 -0.7 0.456
center465 -0.2128 0.8083 0.1679 -1.3 0.205
center479 -0.2484 0.7800 0.2233 -1.1 0.266
center484 0.6843 1.9823 0.3235 2.1 0.034
center486 -0.4359 0.6467 0.2719 -1.6 0.109
center487 -0.0086 0.9915 0.1928 0.0 0.965
center507 -0.5512 0.5763 0.2416 -2.3 0.023
center512 -0.2676 0.7652 0.1979 -1.4 0.176
center520 -0.4155 0.6600 0.2125 -2.0 0.051
center523 -0.0094 0.9906 0.2479 0.0 0.970
center527 -0.0829 0.9204 0.1848 -0.4 0.654
center532 -0.1198 0.8871 0.2115 -0.6 0.571
center534 0.1572 1.1702 0.2389 0.7 0.511
center536 -0.3083 0.7347 0.1904 -1.6 0.105
center588 0.9604 2.6128 0.3354 2.9 0.004
center595 0.0278 1.0282 0.2314 0.1 0.904
center615 0.1760 1.1924 0.2417 0.7 0.466
center620 -0.0371 0.9636 0.1697 -0.2 0.827
center633 0.0690 1.0715 0.2075 0.3 0.739
center640 -0.4037 0.6679 0.1864 -2.2 0.030
center642 -0.5537 0.5748 0.2389 -2.3 0.020
center643 -0.1007 0.9042 0.1613 -0.6 0.532
center645 0.1723 1.1881 0.2076 0.8 0.406
center648 -0.1526 0.8584 0.2968 -0.5 0.607
center656 0.4038 1.4975 0.1736 2.3 0.020
center667 -0.3170 0.7283 0.2075 -1.5 0.127
center670 -0.0583 0.9434 0.3351 -0.2 0.862
center675 0.0308 1.0313 0.1897 0.2 0.871
center688 -0.6747 0.5093 0.2831 -2.4 0.017
center690 -0.5054 0.6032 0.2510 -2.0 0.044
center696 -0.4899 0.6127 0.1945 -2.5 0.012
center700 -0.3198 0.7263 0.1997 -1.6 0.109
center701 -0.6389 0.5279 0.2511 -2.5 0.011
center703 -0.6902 0.5015 0.2115 -3.3 0.001
center708 -0.1915 0.8257 0.2182 -0.9 0.380
center733 -0.1575 0.8543 0.2201 -0.7 0.474
center736 -0.3169 0.7284 0.1916 -1.7 0.098
center742 -0.0582 0.9434 0.2968 -0.2 0.844
center746 0.0610 1.0629 0.1794 0.3 0.734
center749 -0.3705 0.6904 0.2182 -1.7 0.090
center838 -0.5825 0.5585 0.2833 -2.1 0.040
tx_risky_don 0.2629 1.3007 0.0283 9.3 <2e-16
era_tx -0.0522 0.9491 0.0309 -1.7 0.091
tx:three_statusStatus 1B 0.9199 2.5091 0.0445 20.7 <2e-16
tx:three_statusStatus 2 1.3475 3.8477 0.0654 20.6 <2e-16
tx:center7 -0.3549 0.7013 0.2816 -1.3 0.208
tx:center8 0.2563 1.2921 0.3187 0.8 0.421
tx:center18 -0.3720 0.6894 0.3049 -1.2 0.222
tx:center25 0.1620 1.1759 0.2992 0.5 0.588
tx:center33 -0.1296 0.8784 0.2467 -0.5 0.599
tx:center47 -0.4173 0.6588 0.3519 -1.2 0.236
tx:center55 -0.4569 0.6332 0.3531 -1.3 0.196
tx:center63 -0.0981 0.9065 0.4012 -0.2 0.807
tx:center64 -0.5621 0.5700 0.3187 -1.8 0.078
tx:center65 0.2659 1.3046 0.3699 0.7 0.472
tx:center67 -0.6446 0.5249 0.3364 -1.9 0.055
tx:center72 -0.5607 0.5708 0.2645 -2.1 0.034
tx:center78 -0.2863 0.7511 0.2485 -1.2 0.249
tx:center79 -0.5211 0.5939 0.3015 -1.7 0.084
tx:center91 -0.2808 0.7552 0.3045 -0.9 0.357
tx:center94 0.1371 1.1469 0.3446 0.4 0.691
tx:center96 -0.0965 0.9080 0.3266 -0.3 0.768
tx:center108 0.1990 1.2202 0.2972 0.7 0.503
tx:center119 -0.5896 0.5545 0.2824 -2.1 0.037
tx:center128 0.0203 1.0205 0.2929 0.1 0.945
tx:center131 0.1267 1.1351 0.2632 0.5 0.630
tx:center136 0.1143 1.1211 0.2663 0.4 0.668
tx:center141 -0.3803 0.6837 0.2643 -1.4 0.150
tx:center147 0.6130 1.8460 0.4227 1.5 0.147
tx:center149 0.2429 1.2749 0.3588 0.7 0.498
tx:center165 -0.4788 0.6195 0.3702 -1.3 0.196
tx:center171 0.8076 2.2424 0.3213 2.5 0.012
tx:center178 -0.0620 0.9399 0.3007 -0.2 0.837
tx:center183 -0.0253 0.9750 0.2864 -0.1 0.930
tx:center185 -0.6063 0.5453 0.3834 -1.6 0.114
tx:center190 -0.1894 0.8274 0.2774 -0.7 0.495
tx:center195 -0.2354 0.7903 0.3236 -0.7 0.467
tx:center199 0.1528 1.1651 0.3614 0.4 0.672
tx:center201 -0.1768 0.8379 0.2829 -0.6 0.532
tx:center210 -0.0543 0.9472 0.2833 -0.2 0.848
tx:center213 0.1505 1.1624 0.3610 0.4 0.677
tx:center221 -0.5454 0.5796 0.2726 -2.0 0.045
tx:center234 -0.7152 0.4891 0.4663 -1.5 0.125
tx:center298 -0.5603 0.5710 0.2822 -2.0 0.047
tx:center301 -0.2225 0.8005 0.2836 -0.8 0.433
tx:center303 -0.1274 0.8804 0.2875 -0.4 0.658
tx:center307 -0.2421 0.7850 0.3096 -0.8 0.434
tx:center313 -0.5680 0.5667 0.2747 -2.1 0.039
tx:center324 -0.0380 0.9628 0.3177 -0.1 0.905
tx:center337 -0.4626 0.6296 0.2652 -1.7 0.081
tx:center338 -0.7673 0.4643 0.3866 -2.0 0.047
tx:center348 -0.6200 0.5379 0.2823 -2.2 0.028
tx:center350 -0.2140 0.8074 0.2737 -0.8 0.434
tx:center351 -0.4278 0.6519 0.2545 -1.7 0.093
tx:center359 0.1574 1.1705 0.3073 0.5 0.608
tx:center372 -0.4777 0.6202 0.3251 -1.5 0.142
tx:center377 -0.1259 0.8817 0.4304 -0.3 0.770
tx:center380 -0.6783 0.5075 0.3167 -2.1 0.032
tx:center382 -0.0853 0.9183 0.2455 -0.3 0.728
tx:center388 -0.1189 0.8879 0.3119 -0.4 0.703
tx:center402 0.2823 1.3262 0.3266 0.9 0.387
tx:center404 -0.3452 0.7081 0.2658 -1.3 0.194
tx:center408 0.1537 1.1661 0.3451 0.4 0.656
tx:center434 -0.1631 0.8495 0.2387 -0.7 0.494
tx:center442 -0.5739 0.5633 0.3111 -1.8 0.065
tx:center445 -0.4197 0.6573 0.3190 -1.3 0.188
tx:center446 0.6026 1.8269 0.2761 2.2 0.029
tx:center459 -0.1079 0.8978 0.3199 -0.3 0.736
tx:center465 -0.4640 0.6288 0.2468 -1.9 0.060
tx:center479 -0.2044 0.8151 0.3142 -0.7 0.515
tx:center484 -0.6068 0.5451 0.4789 -1.3 0.205
tx:center486 0.7214 2.0573 0.4187 1.7 0.085
tx:center487 -0.1023 0.9027 0.2896 -0.4 0.724
tx:center507 -0.2883 0.7495 0.3294 -0.9 0.382
tx:center512 -0.1300 0.8781 0.2867 -0.5 0.650
tx:center520 0.4900 1.6323 0.2969 1.7 0.099
tx:center523 -0.1633 0.8493 0.3528 -0.5 0.643
tx:center527 -0.2133 0.8079 0.2573 -0.8 0.407
tx:center532 0.0027 1.0027 0.3025 0.0 0.993
tx:center534 -0.3708 0.6902 0.3536 -1.0 0.294
tx:center536 -0.1718 0.8421 0.2602 -0.7 0.509
tx:center588 -1.4317 0.2389 0.4340 -3.3 0.001
tx:center595 -0.6684 0.5125 0.3324 -2.0 0.044
tx:center615 -0.1448 0.8652 0.3329 -0.4 0.663
tx:center620 -0.4204 0.6568 0.2646 -1.6 0.112
tx:center633 0.1392 1.1494 0.2922 0.5 0.634
tx:center640 -0.2954 0.7442 0.2831 -1.0 0.297
tx:center642 0.2081 1.2313 0.3538 0.6 0.556
tx:center643 -0.0308 0.9696 0.2424 -0.1 0.899
tx:center645 -0.5774 0.5614 0.3015 -1.9 0.056
tx:center648 0.4840 1.6225 0.3834 1.3 0.207
tx:center656 -0.6646 0.5145 0.2616 -2.5 0.011
tx:center667 -0.0556 0.9459 0.2921 -0.2 0.849
tx:center670 -0.2930 0.7460 0.4865 -0.6 0.547
tx:center675 -0.2036 0.8158 0.2638 -0.8 0.440
[ reached getOption("max.print") -- omitted 13 rows ]
Likelihood ratio test=3849 on 213 df, p=<2e-16
n= 162054, number of events= 10962
sum_fe_model <- summary(full_fe_model)$coefficients
tx_1a <- sum_fe_model[1,1]
fe_coef <- tibble(name = names(full_fe_model$coefficients),
beta = unname(sum_fe_model[,1])) %>%
filter(grepl( "center",name)== TRUE) %>%
mutate(tx = ifelse(grepl( "tx:",name), "transplant", "waitlist"),
center = numextract(name)) %>%
select(-name) %>%
spread(tx, beta)
vcov_full_fe <- vcov(full_fe_model)
tx_wait_cov <- vector("double", length = length(unique(for_fe$center)))
center_ids <- vector("double", length = length(unique(for_fe$center)))
wait_variance <- vector("double", length = length(unique(for_fe$center)))
tx_variance <- vector("double", length = length(unique(for_fe$center)))
i <- 1
j <- 1
for (n in names(vcov_full_fe[,1])){
if (grepl("center", n) == TRUE & grepl("tx", n) == FALSE){
tx_effect <- paste0("tx:",n)
tx_wait_cov[[i]] <- vcov_full_fe[[tx_effect,n]]
wait_variance[[i]] <- vcov_full_fe[[n,n]]
center_ids[[i]] <- numextract(n)
i <- i + 1
}
if (grepl("tx:center", n) == TRUE){
tx_variance[[j]] <- vcov_full_fe[[n, n]]
j <- j + 1
}
}
tx_wait_cov <- tibble(center = center_ids,
wait_variance = wait_variance,
tx_variance = tx_variance,
covariance = tx_wait_cov)
fe_coef <- fe_coef %>%
left_join(tx_wait_cov) %>%
mutate(tx_benefit = transplant - waitlist) %>%
arrange(-tx_benefit) %>%
mutate(id = row_number())
Joining, by = "center"
ggplot(fe_coef, aes(x = tx_benefit)) + geom_histogram(breaks = seq(-3, 2, 0.25))
labs(x = "center effect on benefit of transplant\n(log hazard ratio of transplant (Status 1A)", y = "number of centers") +
theme_few()
NULL
ggsave("eFigure6.pdf")
s_w_norm_test <- comma_2(shapiro.test(fe_coef$tx_benefit)[[2]])
Distribution of center-specific survival benefits for a Status 1A recipient, calculated by subtracting the log hazard of waitlist from the the log hazard of transplant at the center. The shapiro-wilk test for non-normality was not significant (p = 0.18). The variance of the survival benefit of transplant was 0.378 on log hazard ratio scale, -135% lower than the three-status model.
\(benefit = \beta_{tx, i} - \beta_{wait, i}\)
for_spear <- center_re %>%
filter(center %in% fe_coef$center) %>%
left_join(fe_coef %>% select(center, tx_benefit))
spearman_cor <- cor(for_spear$benefit, for_spear$tx_benefit,method = "spearman")
ggplot(for_spear, aes(x = benefit, y = tx_benefit)) + geom_point() +
labs( x = "benefit of transplant (RE estimate) log hazard",
y = "benefit of transplant (FE estimate) log hazard")
ggsave("eFig7.pdf")
Scatter plot of center fixed-effect (y-axis) and random-effect (x-axis) estimates of the survival benefit of heart transplantation (on log hazard scale). A spearman correlation between the two was -0.9696986
ggplot(fe_coef, aes(x = waitlist, y = tx_benefit)) + geom_point(size = 2.5, alpha = 0.9) +
labs(x = "Waitlist relative risk (log hazard)",
y = "Relative benefit of transplant (log hazard)") +
geom_line(aes(linetype = "Mean Center"),
y = 0, color = "black") +
scale_linetype_manual(values = "dashed")
slope_fe <- lm(tx_benefit ~ waitlist, data = fe_coef)
slope_pct_fe <- comma_2(slope_fe$coefficients[[2]])
ci_slope_fe_wait_low <- comma_2(10*-confint(slope_fe, 'waitlist', level=0.95)[1,1])
ci_slope_fe_wait_up <- comma_2(10*-confint(slope_fe, 'waitlist', level=0.95)[1,2])
ggsave("eFigure8.pdf")
The association between estimated log hazard of death on waitlist at a given center (x-axis) and log hazard survival benefit (y-axis). A linear association is observed, for every one unit increase in the log hazard of waitlist risk the log hazard of transplant decreased -1.8 (95% CI 19 - 16).
to_model <- read.csv("clean_time_series_revise2006_2015.csv") %>%
mutate(
ischemia = case_when(
REC_HR_ISCH/60 < 2 ~ 1,
REC_HR_ISCH/60 >= 2 & REC_HR_ISCH/60 < 4 ~ 2,
REC_HR_ISCH/60 >=4 & REC_HR_ISCH/60 < 6 ~ 3,
REC_HR_ISCH/60 >= 6 & REC_HR_ISCH/60 < 8 ~4,
TRUE ~ 5
),
don_age = factor(case_when(
DON_AGE >=40 & DON_AGE <50 ~ "30-50",
DON_AGE >= 50 ~ "> 50",
TRUE ~ "< 30"
)),
don_renal_function =
case_when(
(DON_BUN/DON_CREAT) > 30 ~ 1,
TRUE ~ 0
),
race_mismatch =
case_when(
black != black_don ~ 1,
TRUE ~ 0
)
)
better_donor_year <- "Surv(t_1, t_2, dead) ~ tx*(stat_2 + stat_1b + factor(list_year)) + don_renal_function + don_age + race_mismatch + ischemia"
six-status, donor quality and transplant year better_donor_year_candidate <- "Surv(t_1, t_2, dead) ~ tx*(stat_2 + stat_1b + factor(list_year) + age + female + bmi_low + bmi_high + simple_diag + cross_match +blood_type) + don_renal_function + don_age + race_mismatch + ischemia"
three status, region & LVAD region_lvad <- "Surv(t_1, t_2, dead) ~ tx*(stat_2 + stat_1b + factor(REGION) + cf_lvad) + tx_risky_don + era_tx"
#three-status, donor quality, and transplant year
better_donor_year <- "Surv(t_1, t_2, dead) ~ tx*(stat_2 + stat_1b + factor(list_year)) + don_renal_function + don_age + race_mismatch + ischemia"
#six-status, donor quality and transplant year
better_donor_year_candidate <- "Surv(t_1, t_2, dead) ~ tx*(stat_2 + stat_1b + factor(list_year) + age + female + bmi_low + bmi_high + simple_diag + cross_match +blood_type) + don_renal_function + don_age + race_mismatch + ischemia"
#three status, region & LVAD
region_lvad <- "Surv(t_1, t_2, dead) ~ tx*(stat_2 + stat_1b + factor(REGION) + cf_lvad) + tx_risky_don + era_tx"
me_models_alt <- vector(mode = "list", length(formulas_alt))
for (i in seq_along(formulas_alt)) {
f_coxme <- as.formula(paste(formulas_alt[[i]], " + (1+ tx|center)"))
start_time <- Sys.time()
fit <- coxme(f_coxme, data = to_model)
me_models_alt[[i]] <- fit
end_time <- Sys.time()
diff_time <- end_time - start_time
print(f_coxme)
print(diff_time)
}
save.image("end_results_with_alts.RData")
print(me_models_alt[[1]])
Cox mixed-effects model fit by maximum likelihood
Data: alt_to_model
events, n = 11052, 163109
Iterations= 24 148
NULL Integrated Fitted
Log-likelihood -107534.9 -105726 -105540.2
Chisq df p AIC BIC
Integrated loglik 3617.96 31.00 0 3555.96 3329.34
Penalized loglik 3989.47 144.03 0 3701.42 2648.55
Model: f_coxme
Fixed coefficients
coef exp(coef) se(coef) z p
tx -2.40330252 0.09041885 0.07950817 -30.23 0.0e+00
stat_2 -1.43060646 0.23916383 0.03726867 -38.39 0.0e+00
stat_1b -0.94614340 0.38823541 0.03166912 -29.88 0.0e+00
factor(list_year)2007 -0.03922706 0.96153236 0.06562071 -0.60 5.5e-01
factor(list_year)2008 -0.10269913 0.90239844 0.06336356 -1.62 1.1e-01
factor(list_year)2009 -0.15341267 0.85777567 0.06214082 -2.47 1.4e-02
factor(list_year)2010 -0.23099307 0.79374497 0.06254102 -3.69 2.2e-04
factor(list_year)2011 -0.31173899 0.73217260 0.06496752 -4.80 1.6e-06
factor(list_year)2012 -0.30452275 0.73747525 0.06386275 -4.77 1.9e-06
factor(list_year)2013 -0.40327413 0.66812892 0.06384416 -6.32 2.7e-10
factor(list_year)2014 -0.33509032 0.71527348 0.06241014 -5.37 7.9e-08
factor(list_year)2015 -0.48253235 0.61721840 0.06590417 -7.32 2.4e-13
don_renal_function -0.03362085 0.96693805 0.02783947 -1.21 2.3e-01
don_age> 50 0.38436954 1.46868808 0.04551777 8.44 0.0e+00
don_age30-50 0.20195758 1.22379609 0.03462113 5.83 5.4e-09
race_mismatch 0.17018309 1.18552189 0.03017527 5.64 1.7e-08
ischemia 0.14592292 1.15710699 0.02361290 6.18 6.4e-10
tx:stat_2 1.40497859 4.07543950 0.06474647 21.70 0.0e+00
tx:stat_1b 0.91921382 2.50731840 0.04384848 20.96 0.0e+00
tx:factor(list_year)2007 -0.03132753 0.96915809 0.08357561 -0.37 7.1e-01
tx:factor(list_year)2008 0.05827196 1.06000323 0.08281571 0.70 4.8e-01
tx:factor(list_year)2009 0.11930700 1.12671577 0.08298896 1.44 1.5e-01
tx:factor(list_year)2010 0.11326307 1.11992651 0.08522732 1.33 1.8e-01
tx:factor(list_year)2011 0.25650368 1.29240352 0.08865318 2.89 3.8e-03
tx:factor(list_year)2012 0.23540186 1.26541719 0.08958282 2.63 8.6e-03
tx:factor(list_year)2013 0.45418073 1.57488260 0.09051088 5.02 5.2e-07
tx:factor(list_year)2014 0.44445880 1.55964589 0.09224301 4.82 1.4e-06
tx:factor(list_year)2015 0.70528065 2.02441477 0.09787104 7.21 5.8e-13
Random effects
Group Variable Std Dev Variance Corr
center Intercept 0.21296954 0.04535602 -0.48052397
tx 0.23821986 0.05674870
write_csv(table_fe_results(me_models_alt[[1]]), "expanded_donor_time.csv")
three_stat_don_year <- variance_survival_benefit(me_models_alt[[1]])
The variance of the survival benefit of transplant was 0.151 on log hazard ratio scale, 6.3% lower than the three-status model.
print(me_models_alt[[2]], digits = 2)
Cox mixed-effects model fit by maximum likelihood
Data: alt_to_model
events, n = 11052, 163109
Iterations= 20 124
NULL Integrated Fitted
Log-likelihood -107534.9 -105459.5 -105275.3
Chisq df p AIC BIC
Integrated loglik 4151 53 0 4045 3657
Penalized loglik 4519 165 0 4189 2979
Model: f_coxme
Fixed coefficients
coef exp(coef) se(coef) z p
tx -1.49188 0.22 0.1260 -11.84 0.0e+00
stat_2 -1.45977 0.23 0.0377 -38.72 0.0e+00
stat_1b -0.93946 0.39 0.0318 -29.57 0.0e+00
factor(list_year)2007 -0.03458 0.97 0.0658 -0.53 6.0e-01
factor(list_year)2008 -0.08925 0.91 0.0636 -1.40 1.6e-01
factor(list_year)2009 -0.13473 0.87 0.0624 -2.16 3.1e-02
factor(list_year)2010 -0.22585 0.80 0.0628 -3.59 3.2e-04
factor(list_year)2011 -0.29964 0.74 0.0653 -4.59 4.5e-06
factor(list_year)2012 -0.28734 0.75 0.0643 -4.47 7.9e-06
factor(list_year)2013 -0.37461 0.69 0.0644 -5.82 5.9e-09
factor(list_year)2014 -0.28896 0.75 0.0632 -4.57 4.8e-06
factor(list_year)2015 -0.45012 0.64 0.0666 -6.76 1.3e-11
age 0.01486 1.01 0.0012 11.94 0.0e+00
female 0.01608 1.02 0.0328 0.49 6.2e-01
bmi_low 0.06362 1.07 0.0360 1.77 7.7e-02
bmi_high 0.05618 1.06 0.0319 1.76 7.8e-02
simple_diagIschemic cardiomyopathy 0.13350 1.14 0.0328 4.07 4.7e-05
simple_diagOther 0.53039 1.70 0.0425 12.48 0.0e+00
simple_diagRestrictive cardiomyopathy 0.19077 1.21 0.0519 3.68 2.3e-04
cross_match 0.39633 1.49 0.0350 11.34 0.0e+00
blood_typeAB 0.13814 1.15 0.0834 1.66 9.7e-02
blood_typeB -0.00089 1.00 0.0461 -0.02 9.8e-01
blood_typeO 0.00624 1.01 0.0300 0.21 8.4e-01
don_renal_function -0.01947 0.98 0.0281 -0.69 4.9e-01
don_age> 50 0.37988 1.46 0.0459 8.27 1.1e-16
don_age30-50 0.19370 1.21 0.0348 5.57 2.5e-08
race_mismatch 0.18680 1.21 0.0307 6.08 1.2e-09
ischemia 0.13947 1.15 0.0236 5.90 3.6e-09
tx:stat_2 1.43691 4.21 0.0665 21.62 0.0e+00
tx:stat_1b 0.91166 2.49 0.0441 20.66 0.0e+00
tx:factor(list_year)2007 -0.03424 0.97 0.0837 -0.41 6.8e-01
tx:factor(list_year)2008 0.04992 1.05 0.0830 0.60 5.5e-01
tx:factor(list_year)2009 0.10972 1.12 0.0832 1.32 1.9e-01
tx:factor(list_year)2010 0.11815 1.13 0.0855 1.38 1.7e-01
tx:factor(list_year)2011 0.25129 1.29 0.0890 2.82 4.8e-03
tx:factor(list_year)2012 0.23022 1.26 0.0900 2.56 1.1e-02
tx:factor(list_year)2013 0.43592 1.55 0.0911 4.78 1.7e-06
tx:factor(list_year)2014 0.41142 1.51 0.0930 4.42 9.7e-06
tx:factor(list_year)2015 0.69276 2.00 0.0986 7.03 2.1e-12
tx:age -0.01629 0.98 0.0018 -9.25 0.0e+00
tx:female 0.00115 1.00 0.0467 0.02 9.8e-01
tx:bmi_low -0.12500 0.88 0.0498 -2.51 1.2e-02
tx:bmi_high 0.12326 1.13 0.0474 2.60 9.4e-03
tx:simple_diagIschemic cardiomyopathy 0.11725 1.12 0.0464 2.53 1.2e-02
tx:simple_diagOther -0.37577 0.69 0.0661 -5.69 1.3e-08
tx:simple_diagRestrictive cardiomyopathy -0.10533 0.90 0.0753 -1.40 1.6e-01
tx:cross_match -0.29432 0.75 0.0585 -5.03 4.9e-07
tx:blood_typeAB -0.13313 0.88 0.1043 -1.28 2.0e-01
tx:blood_typeB 0.06121 1.06 0.0620 0.99 3.2e-01
tx:blood_typeO 0.08479 1.09 0.0432 1.96 5.0e-02
Random effects
Group Variable Std Dev Variance Corr
center Intercept 0.220 0.048 -0.462
tx 0.222 0.049
write_csv(table_fe_results(me_models_alt[[2]]), "expanded_donor_time_can_factors.csv")
three_stat_don_year_can <- variance_survival_benefit(me_models_alt[[2]])
The variance of the survival benefit of transplant was 0.143 on log hazard ratio scale, 11% lower than the three-status model.
print(me_models_alt[[3]])
Cox mixed-effects model fit by maximum likelihood
Data: alt_to_model
events, n = 11052, 163109
Iterations= 38 232
NULL Integrated Fitted
Log-likelihood -107534.9 -105678.7 -105520.1
Chisq df p AIC BIC
Integrated loglik 3712.48 32.00 0 3648.48 3414.55
Penalized loglik 4029.76 107.21 0 3815.33 3031.56
Model: f_coxme
Fixed coefficients
coef exp(coef) se(coef) z p
tx -2.357509614 0.09465566 0.13545664 -17.40 0.000
stat_2 -1.570915146 0.20785488 0.03891259 -40.37 0.000
stat_1b -0.914340893 0.40078070 0.03180455 -28.75 0.000
factor(REGION)2 0.162890733 1.17690809 0.12791626 1.27 0.200
factor(REGION)3 0.102456977 1.10788964 0.12903336 0.79 0.430
factor(REGION)4 0.241820296 1.27356531 0.13127578 1.84 0.065
factor(REGION)5 0.137378133 1.14726188 0.12920832 1.06 0.290
factor(REGION)6 0.110138494 1.11643268 0.17952348 0.61 0.540
factor(REGION)7 0.104057882 1.10966468 0.13217237 0.79 0.430
factor(REGION)8 0.008497164 1.00853337 0.15438771 0.06 0.960
factor(REGION)9 0.033757128 1.03433337 0.15114072 0.22 0.820
factor(REGION)10 0.150832619 1.16280201 0.13754055 1.10 0.270
factor(REGION)11 0.132840629 1.14206797 0.12817764 1.04 0.300
cf_lvad -0.484829417 0.61580223 0.03251422 -14.91 0.000
tx_risky_don 0.264771016 1.30313255 0.02803009 9.45 0.000
era_tx -0.025422473 0.97489796 0.03067013 -0.83 0.410
tx:stat_2 1.632155649 5.11488875 0.06669965 24.47 0.000
tx:stat_1b 0.913485484 2.49299671 0.04427438 20.63 0.000
tx:factor(REGION)2 0.086131081 1.08994919 0.15135504 0.57 0.570
tx:factor(REGION)3 0.126922422 1.13532894 0.15340251 0.83 0.410
tx:factor(REGION)4 0.117464298 1.12464148 0.15528196 0.76 0.450
tx:factor(REGION)5 -0.012114629 0.98795846 0.15270754 -0.08 0.940
tx:factor(REGION)6 -0.515973433 0.59691925 0.21897002 -2.36 0.018
tx:factor(REGION)7 -0.016365998 0.98376720 0.15760212 -0.10 0.920
tx:factor(REGION)8 0.069539339 1.07201423 0.18036698 0.39 0.700
tx:factor(REGION)9 0.071850262 1.07449444 0.17778277 0.40 0.690
tx:factor(REGION)10 0.013153071 1.01323995 0.16333119 0.08 0.940
tx:factor(REGION)11 0.073577763 1.07635223 0.15236926 0.48 0.630
tx:cf_lvad 0.677786619 1.96951362 0.04494613 15.08 0.000
Random effects
Group Variable Std Dev Variance Corr
center Intercept 0.19761376 0.03905120 -0.55829070
tx 0.19753476 0.03901998
write_csv(table_fe_results(me_models_alt[[3]]), "btt_lvad_and_region.csv")
btt_lvad_region <- variance_survival_benefit(me_models_alt[[3]])
The variance of the survival benefit of transplant in the model including UNOS region and LVAD was 0.122 on log hazard ratio scale, 24% lower than the three-status model.